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FORTRAN  PROGRAMS  FOR  COMPUTING  THE  SPACE-TIME 
CORRELATIONS  OF  TURBULENT  FLOW  IN  A  CHANNEL 


1.  Introduction 

I  t. 

r  r  '  / 

Recently,  several  groups  have  performed  direct  numerical  simulations  of  turbulent 

i:  ^ 

flow  in  a  channel  at  low  Reynolds  number s|^  These  have  included  the  high  resolution 
(fully  resolved)  calculation  of  Kim,  Moin,  arid  Moser  (1987)  and  the  low  resoluton  calcu¬ 
lations  of  Gilbert  and  Kleiser  (1987).  'f’hese  simulations  generate  enormous  cjuantities 
of  data  which  have  been  used  to  investigate  the  validity  of  current  theories  of  turbu¬ 
lence,,  to  study  mechanisms  for  the  production  of  turbulence  in  wall-bounded  flows,  and 

C) 

to  develop  new  turbulence  models.  ^However,  before  s^ueh  detailed  investigations  can  be 
undertaken,  it  is  important  to  first  check  some  basic  statistics  of  the  turbulence  such 
as  the  mean  velocity  profile,  turbulence  intensities,  and  Reynolds  stress  profiles.  Be¬ 
yond  these  one-point  statistics,  it  is  also  of  considerable  value  to  have  a  rapid  means  of 
computing  the  two-point  statistics  .so  that  the  length  and  time  scales  of  the  turbulence 
can  be  determined,  dn  this  report,  we  describe  FORTRAN  codes|.which  we  have  devel¬ 
oped  which  enable  rapid  evaluation  of  the  two  point  correlations  from  a  given  dataset.  — ^ 
Datasets  were  crecited  in  two  forms.  In  the  first,  we  have  taken  ‘spatial  snapshots’ 
(realizations)  of  the  flow.  That  is,  at  a  given  time  in  the  calculation,  we  have  stored 
the  three  components  of  the  velocity  at  each  collocation  point  designated  by  its  coor¬ 
dinates  x,y.  and  z.  which  are  respectively  the  streamwise,  spanwise,  and  wall-normal 
c(jordinates.  In  general,  we  have  saved  the  realizations  at  time  intervals  on  the  order 
of  or  larger  than  the  largest  correlation  time  of  the  flow  so  that  we  may  say  that  the 
datasef  <  onsists  of  a  number  of  statistically  independent  snapshots.  In  the  second  type 
Manuscript  approved  October  6,  1988. 
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of  data,  we  have  stored,  at  certain  fixed  locations  above  the  wall,  the  three  components 
of  velocity  in  the  horizontal  planes  at  equally  spaced  time  intervals.  The  time  between 
saves  for  the  planar  data  is  generally  significantly  less  than  the  correlation  time  of  the 
fiow.  We  will  call  such  datasets  ‘planar’  datasets.  In  Section  2  we  describe  codes  which 
produce  the  two-point  correlation  function  of  turbulence  from  a  ‘snapshot’  dateiset  and 
in  Section  3  we  describe  the  code  which  produces  the  space- time  correlation  function 
from  a  ‘planar’  dataset. 

2.  Computation  of  the  Complete  Two-Point  Correlation  Function 
for  Turbulent  Channel  Flow  from  Spatial  Realizations 

A.  Mathematical  Preliminaries  and  Definitions 


Before  describing  the  computer  codes  in  detail,  we  briefly  discuss  certain  mathe¬ 
matical  definitions  and  the  methods  used  to  compute  the  two-point  correlation  function. 
(In  the  following  discussions,  we  use  the  variables  a.’i,X3,  and  X2  interchangeably  with 
the  notation  x,  y,  x,  and  bold  face  notation  indicates  a  vector  quantity.)  For  any  given 
flow  field,  u(x,t),  we  obtain  a  set  of  statistically  independent  realizations  by  sampling 
the  field  at  times  separated  by  at  least  the  characteristic  correlation  time  of  the  flow. 
We  define  the  Fourier  coefficients  of  any  such  realization  by: 


iVi-lNa-l 

$(n,m,.r2)=  ^  u(/, ik, X2)e"^e"^ 

/=0  k=Q 


^  =  n/=M, 
n  =  0,....^  +  1, 


-iV3  iVa  ^ 
m  =  ^— ,....-— -1. 
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The  expression  above  defines  collocation  points  (xi ),  and  and  wavenumbers  (A:i  )„ 

and  (A:3)^  as  follows; 

{xi)i  =  I  X  Li/Ni, 

(^3)*  =  k  X  L3/N3,  (2) 

...  27r 

(Ari)„  =  n  X  — , 

27r 

(^'3)m  =  m  X  — , 

in  which  Ti  and  L3  are  the  domain  lengths  in  the  streamwise  and  spanwise  directions 
respectively.  In  equation  1  time  serves  merely  as  a  parameter  which  may  be  used  to 
identify  any  given  flow  realization  and  is  therefore  suppressed.  We  now  follow  Sirovich 
(1987),  who  suggests  that  the  geometric  symmetries  of  a  flow  should  be  exploited  to 
increase  the  effectix^e  number  of  realizations  that  are  used  in  the  averaging  process.  For 
the  case  of  channel  flow,  the  equations  of  motion  are  invariant  with  respect  to  wall 
normal  reflection,  spanwise  reflection,  and  180  degree  rotation  about  the  streamwise 
axis.  Applying  these  symmetries  to  the  spectral  representation  where  we  use  the 
superscript  0  to  refer  to  the  representation  before  reflections  and  rotations  have  been 
applied,  yields  three  additional  representations  of  the  flow  which  are  given  by: 

=  {$°(n,7?7.,-X2), -$2('b"b-a-’2), 

$2  _  {$0(71, -,77..r2),$2(»b-”7,X2),-$3(n, -777, X2)},  (3) 

$0  =  {$0(77,  -777,  -.7-2),  -777,  -X2),  -<1*3(77,  -777,  -X2)}. 
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In  these  expressions,  the  superscripts  1,  2,  and  3  refer  to  vertical  reflection,  span- 
wise  reflection,  and  streamwise  rotation  symmetries  respectively.  The  subscripts  con¬ 
tinue  to  represent  the  coordinate  directions  as  previously  deflned.  The  (complex)  spec¬ 
trum  'I’a/j  can  now  be  computed  from; 

3  _ 

'Pa^(n,m,X2,X2')  =  (-^$S(n,m,X2)$^(n,m,X2'))>  (4) 

p=0 

where 


aj=  1,2,3. 


The  overbar  designates  complex  conjugate,  and  the  brackets  represent  an  average  over 
all  the  available  realizations.  We  note  that  in  addition  to  the  increase  in  the  effective 
size  of  the  data  base  by  a  factor  of  four,  the  symmetries  have  the  additional  advantage 
of  reducing  the  size  of  also  by  a  factor  of  four.  A  factor  of  two  comes  from  the 
spanwise  symmetry  which  insures  that  is  unique  for  n  >  0  and  m  >  0.  A  second 
factor  of  two  comes  from  the  vertical  reflection  symmetry  which  gives  unique  values  of 
the  spectrum  for  —  1  <  X2  <  0  and  —1  <  X2'  <  1.  The  two  point  correlation  function  , 
can  now  be  obtained  from: 


cll-i  "a._i 


,  ,  /\  T  /  /\  2rmkt 

fc,X2,X2  )  =  A  2^  '^Q0{n,m,X2,X2  )e  e  ^3  , 


(5) 


n=0  N3 

m= - 


where 


K  = 


_ 1 _ 

y/R  OtOc  (0,0,X2)/2/?/3(0,0,X2') 
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We  note  that  care  must  be  taken  in  the  calculation  of  Ra^  since  is  computed  only 
for  n,m  >  0.  That  is,  the  inherent  symmetries  in  must  be  taken  into  account 
when  equation  5  is  used.  In  particular,  the  following  symmetries  must  hold: 

^'ii(n,-m,a:2,a:2)  =  ^ii(n,  m,  X2,  X2), 


^22{n,-m,X2-,x'2)  =  '^22{ri,m,X2,x'2), 

<ifii{n,-m,X2-,x'2)  =  ^i3{n,m,X2-,x'2), 

^'i2(n,-m,a;2,X2)  =  'I'i2(n,  m,  0:2,  X2),  (6) 

'^I3{n,  -m,X2,X2)  =  -^i3(n,?n,X2,a:2), 

'^23(n,-m,X2,X2)  =  -'^23in,m,X2,X2). 

A  detailed  check  on  the  computation  of  was  made  by  computing  the  energies 
(mean  square  values)  of  each  velocity  component  by  summing  over  m  and  n.  The 
comparison  of  these  values  with  the  known  energies  was  exact  and  will  be  described  in 
a  future  report. 

B.  Description  of  the  codes  KLEXP  and  KLANA 

The  computation  of  the  correlation  function  Rq^  defined  by  eqation  5  is  computed 
by  a  two-stage  process.  First,  the  spectral  function  defined  by  equation  4  is 
computed  and  saved.  This  computation  is  performed  using  the  code  KLEXP.  Then 
the  spectral  function  is  used  to  calculate  Ra^  using  the  code  KLANA.  Originally,  the 
interest  was  primarily  to  perform  an  eigenfunction  decomposition  (Lumley(1970))  of 
turbulent  flow  in  a  channel.  To  perform  such  a  decomposition,  the  code  KLEXP  was 
written  to  determine  which  was  of  greatest  interest  for  our  purposes.  Later  it  was 
decided  that  the  the  real  space  function  R^ff  would  also  be  of  interest  so  that  the  code 
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KLANA  was  written  for  the  purpose  of  computing  it.  All  codes  discussed  in  this  report 
are  listed  in  the  Appendix. 

B.l  Description  of  the  code  KLEXP 

The  code  KLEXP  listed  in  the  Appendix  was  written  to  compute  from  a 
data  set  which  consists  of  30  realizations  of  channel  flow  turbulence.  In  this  particular 
case,  there  were  7  data  files  containing  4  realizations  and  one  containing  2.  On  each 
file,  the  three  components  of  vorticity  were  also  stored  but  were  skipped  each  time  the 
data  was  accessed.  After  each  realization  of  the  velocity  field  is  read  into  memory, 
certain  baseline  statistics  such  as  the  mean  and  variance  of  each  velocity  component 
are  computed.  These  values,  which  are  calculated  from  the  raw  (real  space)  data, 
can  later  be  compared  to  the  values  obtained  from  the  spectral  function  as  a 
check  on  the  calculation.  Next,  the  horizontal  transform  of  the  three  components  of 
velocity  is  performed  using  the  subroutine  YXFOUR.  This  calculation  is  represented 
mathematically  by  equation  1.  The  six  components  of  fii'e  then  calculated  from 
the  raw  spectral  representation,  $q,  using  the  subroutine  MATRIX.  In  this  routine,  for 
each  spectral  realization  of  the  flenv,  the  product  is  computed  and  averaging  is 

performed  over  the  symmetry  groups  given  in  equation  3.  We  note  that  the  parameter 
‘is’  is  used  to  distinguish  the  four  possible  ways  (w-hich  will  be  described  below)  in 
which  the  symmetries  need  to  l^e  applied. 

We  now  describe  in  detail  how'  the  six  components  of  'I'q/?  are  computed.  In 
rompuring  4/ , ,  we  first  cfuniJiite  r'|’,  and  ?-},  wdiich  are  given  l)y: 

{ii.in,.r2..r2')  =  ^i(n.  m,  .r2  )^i  (n,  n7,  .r?'),  (7) 


and 


rj ,  ( n.tii.i'z.r/)  =  ^|(?7,  Jn,  —X2  in,  -X2'). 


C 


We  then  apply  the  spanwise  reflection  symmetry  to  rJi  and  rj^  using  the  subroutine 
SYMl.  This  yields  and  which  are  given  by; 

rfi(n,m,X2,X2')  =  rii(n,m,X2,X2')  +  r?j(n,  -m,X2,X2'), 


and 


rii(n,?7j,X2,X2')  =  rli(n,m,X2,X2')  +  r}i(n,  -m,X2,X2').  (8) 

Finally,  the  function  is  computed  from: 

^'ii(n,m,a:2,3;2')  =  rfi(n, 7u, X2, 0:2')  +  (n,  777,  X2,  X2').  (9) 


This  procedure  incorporates  all  three  group  symmetries  given  in  equation  3.  The  same 
procedure  is  used  to  compute  ^^22  and  ^33. 

To  obtain  the  other  three  components  of  the  correlation  matrix,  minor  modifica¬ 
tions  of  the  procedure  outlined  above  are  needed.  In  the  computation  of  ^'12  we  must 
take  into  account  the  change  in  sign  of  the  wall-normal  velocity.  This  is  expreesed  by  : 

^'i2  (  77  ,  777,a;2,X2')  =  rf2(n,m,X2,X2')  -  rj2(n,m,X2,X2').  (10) 


In  the  computation  of  "Pis  we  must  modify  the  relations  in  equation  S  since  there  is  a 
sign  change  in  the  spanwise  velocity  component  when  the  spanwise  reflection  symmetry 
is  used.  The  resulting  expressions  for  Tjj  and  rfj  are  given  by: 


r(;5(77,  ??7,.T2,.r2')  =  (  77  ,  777,  .72 ,  X2 ')  -  rj3(77,  -  777,  .T  ) ,  .72 ' ), 


iiid 


7-13(0,,  777,. 72,. 72')  =  7-}3(77.,  777,.72,.72')  “  7-}3  (77, -77? ,  .73  ,  .72 ' ).  (11) 
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In  calculating  'I'23  must  take  into  account  both  the  sign  changes  in  the  wall-normal 
velocity  component  and  the  spanwise  velocity  component.  This  results  in  the  following 
relations: 

r23(n,772,.r2,.r2')  =  77i, X2, X2')  -  r%in, -m,X2,X2'), 

rl:i(n,m,X2,X2')  =  rl^{n,m,X2,X2')  -  -777, .T2, X2'),  (12) 

and 


'i>2An->n.X2.X2')  =  7-23  (  77. 777, .T2, •r2' )  -  rl;i{n,  m,X2,X2'). 


(13) 


The  final  result  .  j.  is  then  scaled  appropriiitely  and  diagnostics  are  applied  to  check 
the  answer.  A  good  check  is  to  determine  whether  the  energy  in  'I'q/j,  obtained  by 
integrating  over  all  wavcnumbe7's,  is  ecj^ual  to  the  energies  computed  from  the  raw  (real 
space)  fields.  The  final  result  is  then  saved  on  disc.  No  unnecessary  information  has 
been  computed  or  stored  in  this  calculation  and,  as  noted  pi'eviously,  the  application 
of  the  synmietries  reduces  the  storage  requirements  l^y  a  factor  of  four. 


B.2  Description  of  the  code  KL.4NA 

Cofle  KLAXA  uses  the  spectral  results  produced  by  KLEXP  as  described  in 
erpiation  4  to  comptite  the  correlations  7?^,^  described  in  equation  5.  Before  the  main 
calculation,  all  six  spectral  functions  are  read  from  disc.  The  main  loop  (do  8000) 
coinputes  the  correlation  functions  one  at  a  time  and  writes  each  result  to  disc  succe.s- 
sively.  Since  the  calciiltition  of  each  of  the  six  corrrlatioiis  is  identical  except  for  the 
application  of  the  synnnetries  described  in  equation  C.  a  ch'seription  of  the  calculation 
of  /?ii  will  s<'rve  as  ;i  description  for  the  other  five.  For  each  pair  of  indices  ‘iz’  and 
'iz])  which  designate  wall-normal  distanc<'s.  a.  two  diiiK'iisional  array  ‘phi’  is  extracted 
from  the  four  dimensional  spectral  function  d)!!.  Since  the'  resultant  correlation  must 
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be  real,  the  known  symmetry  properties  of  are  applied  to  the  array  ‘phi’  in  subrou¬ 
tine  SETPHI  to  generate  the  array  ‘phisym’.  The  two  dimensional  Fourier  transform  is 
then  performed  in  subroutine  RIJ  to  complete  the  calculation.  Note  that  the  parameter 
•ij'  in  the  call  to  RIJ  determines  which  of  the  symmetries  given  in  equation  6  will  be 
applied.  The  final  result  is  placed  into  the  larger  array  ‘rrij’,  which  is  then  normal¬ 
ized  using  the  proper  K  defined  in  equation  5.  These  correlation  functions  and  the 
normalization  factors  which  have  been  used  are  written  to  a  disk  file  for  later  plotting. 

3.  Calculation  of  the  Space-Time  Correlation  Function  from  Planar  Data 

A.  Mathematical  Definitions 


Before  describing  the  details  of  the  FORTRAN  code  that  has  been  written  to  com¬ 
pute  the  space-time  correlation  function,  we  first  describe  the  method  of  calculation. 
Gi'TH  the  i^lanar  realization  u{xi,X3,t),  we  define  the  following  spectral  representa¬ 
tions: 


nid 


Nx-\Nz-\ 

^  ^  ^  'i  If  Tt  I  t 

1=0  k=0 


2  y  n  It  2nm  ki 
N3 


Nt  -lNt-1 


Oain.j.x^)  =  2^  2_^Ua{l,k,X3)e  e 


2 irn  I  i  2irk]  i 
Nt 


1=0  k=0 


Nq  —  1  Nt  —  1 
^  ^  2 w  mli  2‘irkj  i 

rc»(m,  =  2^  2^  A;,  .iq  )e  ^’3  e 

1=0  k=o 


(14) 


a- =  1,2,3. 


W'c  note  that  X2  is  suppressed  in  equation  14  since  it  is  assumed  to  be  fixed  in  these 
calculations.  We  also  define  the  (discretized)  time  and  frequency  as  follows: 

t,i  =  —n, 

A'< 
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T  is  the  total  time  in  any  given  realization.  In  these  calculations  we  apply  the  spanwise 
symmetry  to  7  and  tt  in  the  following  way: 


7'  =  {7?(n,-m,f),75(n,-m,t),-75(n,-m,t)}, 


and 

=  {7r?(-m,j,t),7r5(-m,j,<),-7r5(-m,j,t)}.  (16) 

where  the  subscripts  and  superscripts  have  been  defined  in  Section  2.  We  can  then 
compute  the  (complex)  spectra  defined  by: 

Va0{n,m)  =  (^  X/ 

‘  k=0  "p=0 
,  A^3-1 _ 

QopinJ)  =  (—  X]  ^ain,m,{x3)k)0^^{n,m,{x3)k)),  (17) 

^  k=o 

2  ^  7^  2  .  ^  _ 

Ua3{mJ)  =  (—  ^ 

■'  *  k=0  ^p=0 

and 

=  1,2,3. 

We  note  again  that  no  symmetry  has  been  applied  in  the  calculation  of  as  indicated 
in  equation  17.  In  summmary,  we  compute  the  appropriate  spectral  functions  by  first 
averaging  over  the  approjjiiate  symmetries,  then  averaging  over  a  coordinate  ,  and 
finally  averaging  over  the  available  realizations. 


We  then  define,  as  before,  the  two  point  correlation  functions  by: 
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AV  =  —  -  — 

vA<^oo(0.0)A‘^/}^(0,0) 

B.  Description  of  code  WA\  SPEC 


(19) 


The  Fortran  code  ll’AVSPEC  permits  the  analysis  of  data  of  the  ‘planar’  type  as 
described  in  the  Introduction.  This  data  consists  of  the  three  components  of  velocity  Ui , 
iio.  and  (/3  ,  the  streamwise.  spauwise,  and  normal  components  of  the  vector  u(.ri ,  xs ,  t). 
The  cofle  listed  in  the  Appendix  was  written  to  compute  correlations  from  a  dataset 
roiisistiiiii;  of  sf'veii  files  in  which  horizontal  planes  of  the  three  components  of  velocity 
and  the  vorticity  are  stored.  Each  file  contains  512  time  steps.  After  each  realization 
I  (i.e.  a  file  consisting  of  512  time  steps)  is  read  from  disc  using  the  subroutine  READ, 

the  computations  performed  in  equations  14-17  are  performed  by  subroutine  HFENS. 

i 
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(Array  ‘tinp’  in  subroutine  READ  is  a  dummy  used  to  skip  over  the  values  of  vorticity 
contained  in  the  dataset.) 

The  input  to  subroutine  HFEXS  is  the  raw  planar  data,  and  the  output  consists 
of  the  updated  (ensemble  averaged)  spectra  given  l)y  T.  0.  and  IT  defined  in  equaton 
17.  We  describe  only  the  calculation  of  F  in  HFENS,  since  the  calculation  of  0,  and 
n  is  virtually  the  same.  The  first  loop  in  HFEXS  (  do  100)  computes  the  average 
over  time  as  indicated  in  equation  17.  The  second  loop  (do  70)  extracts  the  two 
dimensional  array.s  (  the  y-x  array)  from  the  three  dimensional  data  (y-x-t  arrays). 
Each  two  dimensional  array  is  then  brought  into  wavenumber  (A's  —  Aq)  space  using 
the  subroutine  HFFTl.  The  subroutine  ENSA\T  then  computes  the  square  of  the 
magnitude  of  the  Fourier  coefficients  and  applies  the  spanwise  reflection  symmetry. 
This  completes  the  calculation  of  F.  As  indicated  in  equation  17,  no  sjunmetry  is 
applied  in  the  the  calculation  of  0. 

As  a  check  on  the  calculation,  we  observe  that  if  we  sum  F,  0,  and  11  over  wavenum- 
l)er  and  frequency,  the  results  must  l^e  identical.  This  sum  is  the  mean  square  value 
(the  energy)  of  the  velocity  component.  This  serves  as  a  check  on  the  calculations  and 
therefore,  following  the  action  of  HFEXS,  the  total  energy  for  each  of  the  three  spectra 
are  computed  using  the  routines  EXERGl,  ENERG2,  and  ENERG3.  The  results  of 
these  calcu.la^'Vtns.  which  will  be  given  in  a  sub.sequent  report,  show  that  the  energies 
were  equal  to  within  machine  roundoff  error.  After  the  energies  have  been  computed, 
they  are  used  to  normalize  the  sjx'ctral  results. 

The  sijectra  F.  0.  and  H  are  then  used  to  compute  the  space-time  or  space- 
sjjace  correlations  given  in  <’(|uation  IS.  This  is  accomi)lished  by  first  calling  SETPHl. 
SETPH2.  and  SETPH.3.  each  of  which  tirranges  the  spectral  data  in  an  appropriate 
hniii  for  subsefiuent  Fourier  transforintition.  The  Fourier  transforms  are  performed 
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by  the  routines  RIJl,  RIJ2,  and  RIJ3,  the  correlations  are  normalized  as  indicated  in 
equation  18,  and  the  results  are  saved  on  disc. 

Conclusions 

We  have  described  computer  codes  that  have  been  written  to  compute  the  relevant 
correlation  functions  and  their  corresponding  spectral  representations  from  stochastic 
velocity  data  which  has  been  obtained  from  direct  simulations  of  channel  flow  turbu¬ 
lence.  The  codes  we  have  written  can  be  used  to  analyze  data  that  is  in  one  of  two 
forms  :  (1)  spatial  ‘snapshots’  of  the  flow  ,  and  (2)  ’planar’  datasets.  In  the  first  case 
we  compute  the  relevant  correlations  in  a  two  step  process.  First,  we  compute  the 
four-dimensional  spectrum  given  in  equation  4  in  which  we  average  over  all  available 
realizations  and  also  incorporate  the  relevant  symmetries.  This  intermediate  result  is 
of  particular  importance  since  it  forms  the  basis  for  a  Karhunen-Loeve  or  eigenfunction 
expansion  to  which  we  referred  in  Section  2.  This  calculation  is  performed  by  the  code 
KLEXP.  These  results  are  then  used  to  compute  the  (real-space)  correlation  function 
given  in  ec^uation  5  using  the  code  KLANA.  For  the  second  case,  using  data  of  the  ‘pla¬ 
nar’  form,  the  relevant  spectra  given  in  equation  17  and  correlations  given  in  equation 
IS  are  computed  by  implementing  the  code  WAVSPEC. 

As  we  mentioned  in  the  Introduction,  along  with  the  theoretical  interest  that  exists 
in  computing  these  multiclimentional  statistics,  these  statistics  also  serve  as  diagnostic 
tools  which  can  be  used  to  check  the  very  large  data  bases  that  are  generated  from 
direct  numerical  simulations.  A  description  and  some  preliminary  interpretations  of 
the  results  of  these  calculations  will  be  given  in  a  subsequent  report. 
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Appendix 


Here  we  list  the  codes  KLEXP,  KLANA,  and  WAVSPEC  described  in  the  text.  It 
is  hoped  that  the  text  description  together  with  the  annotations  within  each  code  will 
enable  a  user  to  make  modifications  for  his  particular  application.  Note  that  the  fft 
routines  CFFT2,  RCFFT2,  and  CRFFT2  are  library  subroutines  which  are  generally 
available  on  the  Cray-XMP. 
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program  klexp 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C  -  THE  PURPOSE  OF  THIS  CODE  IS  TO  COMPUTE  THE  C 

C  CORRELATION  MATRIX  FOR  A  3-D  TURBULENT  CHANNEL  FLOW.  C 
C  THE  RESULTS  OF  THE  CALCULATION  ARE  USED  AS  INPUT  TO  AN  C 
C  EIGENVALUE  SOLVER  FROM  WHICH  WE  CAN  DETERMINE  THE  KARHUNEN-  C 
C  LOEVE  EXPANSION  FUNCTIONS.  THE  CORRELATION  MATRIX  IS  OF  THE  C 
C  FORM  uiuj (m,n, z, zprime) ,  WHERE  THE  INDICES  I  AND  J  REPRESENT  THE  C 
C  VELOCITY  COMPONENETS  ( 1=STREAMWISE  2=WALL  NORMAL  3=SPANWISE) ,  C 
C  M  IS  THE  SPANWISE  WAVENUMBER,  N,  THE  STREAMWISE,  AND  Z  IS  THE  WALL  C 
C  NORMAL  COORDINATE. NOTE  THAT  THERE  CAN  ONLY  BE  SIX  INDEPENDENT  C 
C  MEMBERS  OF  UiUj ,  AND  ALSO,  BECAUSE  OF  SYMMETRY  CONSIDERATIONS  C 
C  WE  NEED  ONLY  COMPUTE  1/2  OF  EACH  OF  THESE  MEMBERS.  ALSO  NOTE  THAT  C 
C  WE  ONLY  NEED  TO  COMPUTE  FOR  POSITIVE  M  AND  N  .  C 
C  IN  THE  CHANNEL  DATA  USED  HERE  THERE  WERE  16,  64  AND  33  GRID  C 
C  POINTS  IN  THE  1,  2  AND  3  DIRECTIONS  RESPECTIVELY . IN  SPECTRAL  SPACE  C 
C  THERE  ARE  7  NON-ZERO  STREAMWISE  WAVENUMBERS  AND  23  SPANWISE  C 
C  MODES.  THE  TRUNCATION  IS  DUE  TO  THE  USE  OF  DEALIASING  IN  THE  C 
C  ORIGINAL  CHANNEL  FLOW  CALCULATION.  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
parameter (nz=32 , ny=64 , nx=16) 

parameter ( nzl=nz+l , nzml=nz-l , nzd21=nz/2+l , nzd2=nz/2 ) 
parameter (nyd2=ny/2 , nyd21=ny/2+l) 
parameter (nxp2=nx+2 , nxd2=nx/2 ,nxd21=nx/2+l) 
parameter (nfiles=8) 
common/param/fac,  fad 

dimension  ul (nzl, ny , nxp2) ,u2 (nzl, ny , nxp2) , u3 (nzl, ny , nxp2) 
dimension  oml(nzl,l,l) 

dimension  ulul (nyd2,nxp2 ,nzd21,nzl) ,u2u2 (nyd2 , nxp2 , nzd21, nzl) 
dimension  u3u3 (nyd2 , nxp2 , nzd2l, nzl) ,ulu2 (nyd2,nxp2,nzd21,nzl) 
dimension  ulu3 (nyd2 , nxp2 ,nzd21,nzl) ,u2u3 (nyd2 , nxp2 , nzd21 , nzl) 
dimension  ulm(nzl) ,u2m(nzl) ,u3m(nzl) 
dimension  ulrms(nzl) ,u2rms(nzl) ,u3rms(nzl) 
dimension  ul2rms(nzl) ,ul3rms(nzl) ,u23rms(nzl) 
dimension  si (nzd21) ,s2 (nzd21) ,s3 (nzd21) 
dimension  sl2 (nzd21) , sl3 (nzd21) , s23 (nzd21) 

C  SET  CERTAIN  USEFUL  SCALING  PARAMETERS 

call  setup 

nsym  =  4 
nens  =  nsym*30 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  FOR  INPUT  WE  HAVE  7  FILES  EACH  OF  WHICH  CONTAINS  4  INDEPENDENT 
C  REALIZATIONS  OF  THE  FLOW  AND  AN  8TH  FILE  THAT  HAS  TWO  REALIZATIONS. 

C  THE  FIRST  28  REALIZATIONS  ARE  SPACED  100  VISCOUS  TIME  UNITS 
C  APART  AND  THE  LAST  TWO,  ARE  50  VISCOUS  TIME  UNITS  APART. 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

do  1  nrec  =  l,nfiles 

lu  =  10+(nrec-l) 
if  (nrec  .le.  7)  then 
ntimes  =  4 
else 

ntimes  =  2 

endif 

do  1  nread  =  1, ntimes 

read(lu) ( ( (ul(k, j ,i) ,k=l,nzl) , j=l,ny) , i=l,nx) 
read(lu) (((u3(k,j,i) ,k=l,nzl) ,j=l,ny) ,i=l,nx) 
read(lu) (((u2(k,j,i) ,k=l,nzl) ,j=l,ny) ,i=l,nx) 
read(lu) ( ( (oml (k, j , i) ,k=l,nzl) ,j=l,l) ,i=l,l) 
read(lu) ( ( (oml (k, j , i) ,k=l,nzl) , j=l,l) ,i=l,l) 
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read(lu)  ( ( (oinl(k,  j  ,  i)  ,k=l/nzl)  ,i=l,l) 

CCCCCCCCCCCC  COMPUTE  MEAN  VALUES  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCC  THIS  IS  BEING  DONE  MERELY  FOR  DIAGNOSTIC  PURPOSES.  THAT  IS, 

CC  WE  ARE  COMPUTING  THE  MEAN  AND  RMS  FOR  EACH  OF  THE  THREE  FIELDS, 

call  itiean(ul,ulm) 
call  mean(u2  ,u2in) 
call  inean(u3,u3m) 
call  zms(ul,ul,ulrins) 
call  Zins(u2,u2,u2nns) 
call  zms(u3,u3,u3rins) 
call  zms  (ul,u2  , ul2inns) 
call  zins(ul,u3  ,ul3rins) 
call  zms  (u2  ,  u3  ,  u23rTns) 

CCCCCCCCCCCCC  GO  TO  fOURIER  SPACE  CCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCC  TRANSFORM  EACH  OF  THE  FIELDS  FROM  U(Z,Y,X)  TO  U(Z,M,N) 

CCCC  BUT  FIRST  LETS  INITIALIZE  THE  FFTS. 

call  expo 
call  yxfour(ul) 
call  yxfour(u2) 
call  yxfour(u3) 

CCCCCCCCCCCC  COMPUTE  CORRELATION  MATRICRES  CCCCCCCCCCCCCCCC 
C  COMPUTE  EACH  OF  THE  SIX  INDEPENDENT  MEMBERS  OF  THE  CORRELATION 
C  MATRIX 

call  matrix (ul, ul, ulul, 1) 
call  matrix (u2,u2,u2u2,l) 
call  matrix(u3,u3,u3u3,l) 
call  matrix(ul,u2,ulu2,2) 
call  matrix(ul,u3,ulu3,3) 
call  matrix (u2,u3,u2u3, 4) 

1  continue 

CCCCCCCCCCCCCCCC  NOW  SCALE  RESULTS  CCCCCCCCCCCCCCCCCCCCCCCCCCC 
con  =  1./ (float (30) *fac) 
coni  =  l./(float(nens)*facl) 
call  scale (ulul, coni) 
call  scale (u2u2 , coni) 
call  scale (u3u3 , coni) 
call  scale(ulu2 ,conl) 
call  scale(ulu3 , coni) 
call  scale(u2u3 ,conl) 
call  mult (ulm, con) 
call  mult (u2m, con) 
call  mult (u3m, con) 
call  mult (ulrms, con) 
call  mult (u2rms, con) 
call  mult (u3rms, con) 
call  mult (ul2rms, con) 
call  mult (ul3rms, con) 
call  mult (u23rms, con) 

CCCCCCCCCCCCCCCCCC  COMPUTE  RMS  VALUES  CCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
do  2  k  =  l,nzl 

ulrms(k)  =  sqrt (ulrms (k)  -  ulm(k)**2) 
u2rms(k)  =  sqrt (u2rms(k)  -  u2m(k)**2) 
u3rms(k)  =  sqrt (u3rms (k)  -  u3m(k)**2) 
ul2rms(k)  =  (ul2rms(k)-  ulm(k) *u2m(k) ) 
ul3rms(k)  =  (ul3rms(k)-  ulm(k) *u3m(k) ) 
u23rms(k)  =  (u23rms(k)-  u2m(k) *u3m(k) ) 

2  continue 
ul2rms(l)  =  0,0 
ul2rms(nzl)  =  0.0 
ul3rms(l)  =  0.0 
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ul3rms(nzl)  =  0.0 
u23rins(l)  =  0.0 
u23rms(nzl)  =  0.0 
do  3  k  =2,nzl-l 

ul2rms(k)  =  ul2rms (k) /  (ulrms (k)  *u2rins  (k) ) 
ul3rins(k)  =  ul3rnis(k)/(ulrins(k)  *u3rms(k) ) 
u23rns(k)  =  u2  3rins  (k)/ (u2rins  (k)  *u3rms  (k) ) 

3  continue 

CCCCCCCCCCCCC  COMPUTE  RMS  VALUES  FROM  SPECRAL  RESULTS  CCCCCCCCCCCCCCC 
C  HERE  WE  INTEGRATE  THE  CORRELATION  FUNCTIONS  (AFTER  HAVING  SCALED 

C  THEM  PROPERLY)  TO  SEE  IF  WE  PERFORMED  THE  CALCULATION  PROPERELY. 

C  THAT  IS,  DOES  THE  RMS  WE  GET  FROM  THE  SPECTRA  EQUAL  THE  TRUE  RMS 

C  FOR  EACH  FIELD? 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

call  intg (ulul , si) 
call  intg (u2u2 , s2) 
call  intg (u3u3 , S3 ) 
call  intg (ulu2 , sl2) 
do  4  k  =l,nzd21 

sl(k)  =  sqrt(2.0  *  sl(k)) 
s2(k)  =  sqrt(2.0  *  s2(k)) 

4  S3(k)  =  sqrt(2.0  *  s3(k)) 

sl2(l)  =  0.0 
do  5  k  =  2,nzd21 

sl2(k)  =  (2.*  sl2(k) )/(sl(k)*s2(k) ) 

5  continue 

CCCCCCCCCCCCCCCCCC  WRITE  RESULTS  TO  DISC  CCCCCCCCCCCCCCCCCCCCCC 
C  FIRST  DUMP  THE  CORRELATION  MATRICES 

C  NOTE  WE  ARE  NOW  DUMPING  FROM  IZ  =1,NZD21 

C  AND  IZP  =  1,NZ1.  THAT  IS  ONE  HALF  OF  EACH  MATRIX. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
do  10  izp  =l,nzl 

write (50) ( ( (ulul (j , i, iz, izp) , j=l,23) ,i=l,14) ,iz=l,nzd21) 
write (51) ( ( (u2u2 ( j , i, iz,izp) , j=l,23) , i=l, 14) , iz=l,nzd21) 
write (52) ( ( (u3u3 ( j , i, iz, izp) , j=l,23) , i=l, 14) , iz=l,nzd21) 
write (53) ( ( (ulu2 ( j , i, iz, izp) , j=l, 23) , i=l, 14) , iz=l,nzd21) 
write(54) ( ( (ulu3 ( j , i, iz, izp) ,j=l,23) ,i=l,14) ,iz=l,nzd21) 
write (55) ( ( (u2u3 (j , i, iz, izp) , j=l,23) , i=l, 14) , iz=l,nzd21) 

10  continue 

C  THEN  DUMP  THE  MEAN  VALUES  AND  THE  RMS  VALUES  OF  EACH  FIELD 

write(60, 100)  (ulin(k)  ,k=l,nzl) 
write(61, 100)  (u2in(k)  ,k=l,nzl) 
write(62 , 100)  (u3in(k)  ,k=l,nzl) 
write(63,100) (ulrms(k) ,k=l,nzl) 
write (64 , 100)  (u2rins(k)  ,k=l,nzl) 
write(65, 100)  (u3rins(k)  ,k=l,nzl) 
write (66 , 100) (ul2rms(k) ,k=l,nzl) 
write (67, 100)  {ul3rins(k)  ,k=l,nzl) 
write(68, 100)  (u23rins(k)  ,k=i,nzl) 

CCC  THEN  DUMP  THE  MEAN  SQARE  VALUES  COMPUTED  FROM  THE  SPECTRA. 

C  NOTE  THAT  SINCE  THE  SYMMETRIES  WERE  USED  TO  IN  COMPUTING 
C  THESE  SPECTRA  THE  DEGENERACIES  MUST  BE  INCLUDED  IN  THE 
C  COMPUTATION  OF  THE  ENERGIES.  THIS  HAS  NOT  BEEN  DONE  HERE 

C  SO  THAT  THE  ENERGY  VALUES  WE  GET  HERE  WILL  NOT  AGREE 

C  WITH  THE  VALUES  WE  GET  FROM  THE  RAW  DATA.  THESE  DEGENERACIES 

C  CAN  BE  ACCOUNTED  FOR  AND  THE  RESULTS  WHEN  WE  DO  INCLUDE  THEM 

C  PROPERLY  DO  APEAR  TO  BE  CORRECT. 

write (70, 100) (sl(k) ,k=l,nzd21) 
write (71, 100) (s2 (k) ,k=l,nzd21) 
write(72, 100) (s3(k) ,k=l,nzd21) 
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write (73 ,100) (sl2 (k) ,k=l,nzd21) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
100  format (2x, 4el2 . 5) 
stop 
end 

subroutine  matrix(a,b,c, is) 
parameter (nz=32 , ny=64 , nx=16) 

parameter (nzl=nz+l , nzml=nz-l , nzd2 l=nz/2+l , nzd2=nz/2 ) 
parameter (nyd2=ny/2 , nyd21=ny/2+l) 
parameter (nxp2=nx+2 , nxd2=nx/2 ,nxd21=nx/2+l) 

dimension  a(nzl,ny,nxp2) ,b(nzl,ny,nxp2) , c (nyd2 , nxp2 , nzd21 , nzl) 
dimension  rl(ny,nxp2) ,r2(ny,nxp2) 

C  IN  THIS  SUBROUTINE  WE  COMPUTE  THE  CORRELATION  (C)  FROM  THE 

C  FIELDS  A  AND  B.  WE  UTILIZE  ALL  THE  AVAILABLE  SYMMETRIES  (4)  FOR 

C  THE  CHANNEL  FLOW  PROBLEM.  FOR  EACH  PAIR  OF  WALL  NORMAL  INDICES 

C  IZ  AND  IZP  WE  COMPUTE  THEIR  COUNTERPARTS  IZl  AND  IZ2 .  THE  INDEX 

C  IS  DETERMINES  THE  SIGNS  TO  BE  USED  WHEN  THE  SYMMETRIES  ARE  APPLIED. 

C  THE  SUBROUTINES  SYMl  AND  SYM2  ARE  USED  TO  APPLY  THE  REFLECTIONAL 

C  SYMMETRIES.  PLEASE  NOTE!!!  IZP  GOES  FROM  1  TO  NZl. 

do  1  iz  =  l,nzd21 

izl  =  nzl  -  iz  +  1 
do  1  izp  =  l,nzl 

iz2  =  nzl  -  izp  +  1 

do  2  j=l,ny 
do  2  i=l,nxd21 

ir  =  2*i  -1 
im  =  2*i 

wlr  =  a(iz,j,ir) 
wlim  =  a(iz, j ,im) 
w2r  =  b(izp,j,ir) 
w2im  =  b(izp,j,im) 
wlpr  =  a(izl, j ,ir) 
wipim  =  a(izl,j,im) 
w2pr  =  b ( iz2 , j , ir) 
w2pim  =  b(iz2,j,im) 
rl(j,ir)  =  wlr*w2r  +  wlim*w2im 
rl(j,im)  =  wlr*w2im  -  wlim*w2r 
r2(j,ir)  =  wlpr*w2pr  +  wlpim*w2pim 
r2(j,im)  =  wlpr*w2pim  -  wlpim*w2pr 

2  continue 

CCCCCCCCCCCCCCCCCCCC  INVOKE  SYMMETRIES  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
if  (  is  .eq.  1  )  then 

call  syml (rl,ny,nxp2) 
call  syml (r2 ,ny,nxp2) 
do  3  i=l,nxp2 
do  3  j=l,nyd2 

3  rl(j,i)  =  rl(j,i)  +  r2(j,i) 

endif 

if  (  is  .eq.  2  )  then 

call  syml (rl,ny,nxp2) 
call  syml(r2,ny,nxp2) 
do  4  i=l,nxp2 
do  4  j=l,nyd2 

4  rl(j,i)  =  rl(j,i)  -  r2(j,i) 

endif 

if  (  is  .eq.  3  )  then 

call  sym2 (rl,ny,nxp2) 
call  sym2 (r2,ny,nxp2) 
do  5  i=l,nxp2 
do  5  j=l,nyd2 
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5  rl(j,i)  =  rl(j,i)  +  r2(j,i) 

end  if 

if  (  is  .eq.  4  )  then 

call  syin2  (rl,ny,nxp2) 
call  syin2  (r2,ny,nxp2) 
do  6  i=l,nxp2 
do  6  j=l,nyd2 

6  rl(j,i)  =  rl(j,i)  -  r2(j,i) 

end  if 

CCCCCCCCCCC  CALCULATE  ENSEMBLE  AVERAGE  CCCCCCCCCCCCCCCCCCCCC 
do  7  i  =  l,nxp2 
do  7  j  =  l,nyd2 

c(j , i, iz, izp)  =  c(j,i,iz,izp)  +  rl(j,i) 

7  continue 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

1  continue 

return 
end 

subroutine  setup 

parameter ( ny=64 , nx=16 ) 

parameter (nxp2=nx+2 , nxd21=nx/2+l) 

parameter ( nyt2=2*ny , nyd2 l=ny/2+l , nyd22=ny/2+2 ) 

common/param/fac,  fad 


cccccc  -  some  parameters 

fac  =  float (nx*ny) 
fad  =  ((fac**2)*2) 


return 

end 

subroutine  hfft(a,is) 
parameter (ny=6 4 , nx=16) 
parameter (nxp2=nx+2 , nxd21=nx/2+l) 
parameter ( nyt2=2  *ny , nyd2 l=ny/2+l , nyd22=ny/2+2 ) 
parameter (nx3p4=3*nx+4 , ny5=5*ny , nxny= (nx+2) *ny) 
dimension  a(ny,nxp2) 
common/ index/  indx ( nxp2 ) , indy ( ny t 2 ) 
common/exp/  exrc(nx3p4) ,excr(nx3p4) ,eycc(ny5) 
common/work/  vecx(nxp2) ,vecy(nyt2) ,vecyl(nyt2) 
if  (is.eq.  -1)  go  to  500 
do  1  j=l,ny 
do  2  i=l,nxp2 

indx ( i) = ( i-1) *ny+j 
2  continue 

call  gather (nx,vecx, a, indx) 
call  rcfft2 (0, l,nx,vecx,exrc,vecx) 
call  scatter (nxp2 , a , indx, vecx) 

1  continue 

do  3  i=l,nxd21 
do  4  j=l,ny 

jl=2*(j-l)+l 
indy ( j 1) =j+2* (i-1) *ny 
indy ( j 1+1) =indy ( j 1) +ny 
4  continue 

call  gather (nyt2 ,vecy, a, indy) 
call  cfft2 (0, l,ny,vecy,eycc,vecy) 
call  scatter (nyt2 , a , indy ,vecy) 

3  continue 
return 

500  continue 

do  5  i=l,nxd21 
do  6  j=l,ny 
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jl=2*(j-l)+l 
indy (jl)=j+2* (i-1)  *ny 
indy ( j 1+1 ) =indy ( j 1 ) +ny 
continue 

call  gather (nyt2,vecy, a, indy) 
call  cfft2 (0 ,-l,ny ,vecy ,eycc,vecy) 
call  scatter (nyt2, a, indy, vecy) 

5  continue 

do  7  j=l,ny 
do  8  i=l,nxp2 

indx ( i) = ( i-1) *ny+j 
8  continue 

call  gather (nxp2,vecx, a, indx) 
call  crfft2 (0,-l,nx,vecx,excr,vecx) 
call  scatter (nx, a, indx, vecx) 

7  continue 

do  10  i=l,nx 
do  10  j=l,ny 

a(j , i)=a(j , i)/(2.*float(nx*ny) ) 

10  continue 

return 
end 

subroutine  syml (a , ny , nx) 
dimension  a(ny,nx) 
nyd2=ny/2 
nyd21  =ny/2+l 
do  1  i=l,nx 
do  1  j=l,nyd2 

if  (  j  .eq.  1)  then 

a(j,i)  =  2.0*a{j,i) 
else 

a(j,i)  =  (a(j,i)  +  a(ny-j+2,i)) 
endif 

1  continue 

return 
end 

subroutine  sym2 (a, ny, nx) 

CCCCCC  NOTE:  BECAUSE  OF  THE  SYMMETRY  IN  THIS  CASE  THE  ZERO  SPANWISE 
CCCCCC  MODES  MUST  BE  ZERO 
dimension  a(ny,nx) 
nyd2=ny/2 
nyd21  =ny/2+l 
do  1  i=l,nx 
do  1  j=l,nyd2 

if  (  j  .eq.  1)  then 
a(j,i)  =  0.0 
else 

a(j,i)  =  (a(j,i)  -  a(ny-j+2,i)) 
endif 

1  continue 

return 
end 

subroutine  scale(a,con) 
parameter (nz=3 2 , ny=64 , nx=16) 

parameter (nzl=nz+l , nzml=nz-l , nzd21=nz/2+l , nzd2=nz/2 ) 
parameter ( nyd2=ny/2 , nyd2 l=ny/2+l) 
parameter (nxp2=nx+ 2 , nxd2=nx/2 , nxd21=nx/2+l) 
dimension  a (nyd2 , nxp2 , nzd21 , nzl) 
do  1  izp  =  l,nzl 
do  1  iz  =  l,nzd21 
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a(l, 1, iz, izp)  =0.0 
do  1  i=l,nxp2 
do  1  j=l,nyd2 

1  a ( j , i, iz, izp)  =  a( j , i, iz, izp) *con 

return 
end 

subroutine  mult (a, con) 
parameter (nz=32 , ny=64 , nx=16) 

parameter (nzl=nz+l , nzml=nz-l , nzd2l=nz/2+l , nzd2=nz/2 ) 
parameter (nyd2=ny/2 , nyd21=ny/2+l) 
parameter ( nxp2=nx+2 , nxd2=nx/2 , nxd2 l=nx/2+l ) 
dimension  a(nzl) 
do  1  iz  =  l,nzl 

1  a(iz)  =  a(iz)*con 
return 

end 

subroutine  yxfour(a) 
parameter (nz=32 , ny=64 , nx=16) 

parameter ( nz l=nz+l , nzml=nz-l , nzd2 l=nz/2+l , nzd2=nz/2 ) 
parameter ( nyd2=ny/ 2 , nyd21=ny/2+l) 
parameter (nxp2=nx+2 , nxd2=nx/2 , nxd21=nx/2+l) 
dimension  a (nzl , ny , nxp2) ,work(ny,nxp2) 
do  1  k  =l,nzl 

do  2  j  =l,ny 
do  2  i  =l,nx 

2  work(j,i)  =  a(k, j,i) 
call  hfft(work,l) 

do  3  j  =l,ny 
do  3  i  =l,nxp2 

3  a(k,j,i)  =  work(j,i) 

1  continue 

return 

end 

subroutine  mean(a,s) 
parameter (nz=32 , ny=64 , nx=16) 

parameter (nzl-nz+1, nzml=nz-l, nzd21=nz/ 2+1 , nzd2=nz/2) 
parameter (nyd2=ny/2 , nyd21=ny/2+l) 
parameter (nxp2=nx+2 , nxd2=nx/2 , nxd21=nx/2+l) 
dimension  a (nzl , ny , nxp2) , s (nzl) 
do  2  i  =l,nx 
do  2  j  =l,ny 
do  2  k  =l,nzl 

2  s(k)  =  a(k,j,i)  +  s(k) 
return 

end 

subroutine  zms(a,b,s) 
parameter (nz=3 2 , ny=64 , nx=16) 

parameter (nzl=nz+l , nzml=nz-l, nzd21=nz/2+l ,nzd2=nz/2) 
parameter (nyd2=ny/2 ,nyd21=ny/2+l) 
parameter (nxp2=nx+2 , nxd2=nx/2 , nxd21=nx/2+l) 
dimension  a (nzl , ny , nxp2) ,b(nzl,ny,nxp2) ,s(nzl) 
do  2  i  =l,nx 
do  2  j  =l,ny 
do  2  k  =l,nzl 

2  s(k)  =  a(k, j,i)*b(k, j,i)  +  s(k) 

return 
end 

subroutine  intg(a,s) 
parameter (nz=3 2 ,ny=64 , nx=16) 

parameter (nzl=nz+l, nzml=nz-l,nzd21=nz/2+l,nzd2=nz/2) 
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parameter ( nyd2=ny/2 , nyd2 l=ny/2+l) 
parameter (nxp2=nx+2 , nxd2=nx/2 , nxd21=nx/2+l) 
dimension  a(nyd2,nxp2,nzd21,nzl) ,s(nzd21) 
do  1  iz  =  l,nzd21 
izp  =  iz 
do  1  i=l,nxd21 
do  1  j=l,nyd2 
ir  =  2*i-l 

1  s(iz)  =  a(j ,ir,iz,izp)  +  s(iz) 

return 
end 

subroutine  expo 
parameter (ny=64 , nx=16) 
parameter (nxp2=nx+2 , nxd21=nx/2+l) 
parameter (nyt2=2*ny, nyd21=ny/2+l, nyd22=ny/ 2+2) 
parameter (nx3p4=3  *nx+4 , ny5=5*ny , nxny= (nx+2 ) *ny ) 
common/work/vecx(nxp2) ,vecy(nyt2) ,vecyl(nyt2) 
common/exp/exrc(nx3p4) ,excr(nx3p4) ,eycc(ny5) 
call  rcfft2 (1, l,nx,vecx,exrc, vecx) 
call  crfft2 (1, l,nx,vecx,excr,vecx) 
call  cf ft2 (1, 1 , ny , vecy , eycc , vecy ) 
return 
end 
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program  klana 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C  As  described  in  the  text,  this  code  computes  the 
C  spatial  correlation  function  by  Fourier  transformation  of  the 
C  spectral  correlations  computed  from  the  raw  velocity  fields 
C  by  using  the  klexp  code. 

C 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
parameter ( nz=3  2 , ny=64 , nx=16 ) 
parameter (nym=2 3 , nxm=14) 

parameter ( nz l=nz+l , nzml=nz-l , nzd2 l=nz/2+l , nzd2=nz/2 ) 
parameter (nyd2=ny/2 , nyd21=ny/2+l, nymt2=nym*2) 
parameter (nxp2=nx+2 , nxd2=nx/2 , nxd2 l=nx/2+l , nxmd2=nxm/2 ) 
dimension  z(nzl) ,cnorm(nzd21) ,cn(3,nzl) 

dimension  ulul (nym, nxm, nzd21, nzl) ,u2u2 (nym,nxm,nzd21,nzl) 
dimension  u3u3 (nym,nxm,nzd21,nzl) ,ulu2 (nym,nxm,nzd21,nzl) 
dimension  ulu3 (nym,nxm,nzd21,nzl) ,u2u3 (nym,nxm,nzd21,nzl) 
dimension  phi(nym,nxm) ,phisym(ny,nxp2) ,r(nyd2,nx) 
dimension  rrij (nyd2 , nx,nzd21,nzl) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  First  read  in  the  six  components  of  the  spectral  correlation 
c  function. 

do  10  izp  =  l,nzl 

read (10) ( ( (ulul( j , i, iz, izp) , j=l,nym) , i=l,nxm) , iz=l,nzd21) 
read (11) ( ( (u2u2 ( j , i, iz, izp) , j=l,nym) , i=l,nxm) , iz=l,nzd21) 
read ( 12 ) ( ( (u3u3 ( j , i, iz , izp) , j=l, nym) , i=l, nxm) , iz=l,nzd21) 
read (13) ( ( (ulu2 (j , i, iz, izp) , j=l,nym) , i=l,nxm) , iz=l,nzd21) 
read(14) ( ( (ulu3 ( j , i, iz, izp) , j=l,nym) , i=l,nxm) , iz=l,nzd21) 
read (15) ( ( (u2u3 (j,i,iz,izp) , j=l, nym) , i=l, nxm) , iz=l,nzd21) 

10  continue 

DO  8000  IJ=1,6  I  loop  on  11,22,33,12,13,23  &  calculate  rrij 

DO  5010  IZP=1,NZ1 
do  5010  iz=l,nzd21 

GOTO  (5021,5022,5023,5024,5025,5026),  IJ 

5021  do  5031  j=l,nym 
do  5031  i=l,nxm 

phi(j,i)=ulul(j,i,iz,izp) 

5031  CONTINUE 
goto  5029 

5022  do  5032  j=l,nym 
do  5032  i=l,nxm 

phi ( j , i) =u2u2 ( j , i, iz , izp) 

5032  CONTINUE 
goto  5029 

5023  do  5033  j=l,nym 
do  5033  i=l,nxm 

phi ( j , i) =u3u3 ( j , i , iz, izp) 

5033  CONTINUE 

goto  5029 

5024  do  5034  j=l,nym 
do  5034  i=l,nxm 

phi(j,i)=ulu2(j,i,iz,izp) 
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o  o  o 


5034 

CONTINUE 
goto  5029 

5025 

do  5035  j=l,nym 
do  5035  i=l,nxm 

phi(j,i)=ulu3(j,i,iz,izp) 

5035 

CONTINUE 
goto  5029 

5026 

do  5036  j=l,nym 
do  5036  i=l,nxm 

phi(j , i)=u2u3 (j , i,iz,izp) 

5036 

CONTINUE 

5029 

continue 

call  setphi (phi , phisym, IJ) 
call  rij (phisym, r, cn, ij , iz, izp) 
c 

c  NOW  PUT  r  BACK  INTO  rrij  -  ("rij"  is  name  of  subroutine) 
c 

do  5030  i=l,nx 
do  5030  j=l,nyd2 

rrij(j,i,iz,izp)=r(j,i) 

5030  continue 
5010  continue 
6000  CONTINUE 

Normalize  each  rrij ,  then  write  it  to  disk: 

GOTO  (7501,7501,7501,7504,7505,7506),  IJ 

7501  DO  7511  IZP=1,NZ1 
DO  7511  IZ=1,NZD21 
DO  7511  1=1, NX 
DO  7511  J=1,NYD2 

RRIJ(J,I,IZ,IZP)=RRIJ(J,I,IZ,IZP)/SQRT(CN(IJ,IZ) *CN(IJ,IZP) ) 
7511  CONTINUE 
GOTO  7600 

7504  DO  7514  IZP=1,NZ1 
DO  7514  IZ=1,NZD21 
DO  7514  1=1, NX 

DO  7514  J=1,NYD2 

RRIJ(J,I,IZ,IZP)=RRIJ(J,I,IZ,IZP)/SQRT(CN(1,IZ) *CN(2,IZP) ) 

7514  CONTINUE 
GOTO  7600 

7505  DO  7515  IZP=1,NZ1 
DO  7515  IZ=1,NZD21 
DO  7515  1=1, NX 

DO  7515  J=1,NYD2 

RRIJ(J,I,IZ,IZP)=RRIJ(J,I,IZ,IZP)/SQRT(CN(1,IZ) *CN(3,IZP) ) 

7515  CONTINUE 
GOTO  7600 

7506  DO  7516  IZP=1,NZ1 
DO  7516  IZ=1,NZD21 
DO  7516  1=1, NX 

DO  7516  J=1,NYD2 
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RRIJ(J,I, IZ,IZP)=RRIJ(J,I,IZ,IZP)/SQRT(CN(2,IZ) *CN(3,IZP) ) 
7516  CONTINUE 

7600  CONTINUE 
C 

c  DUMP  rrij  TO  DISKFILE  CONSISTENT  WITH  uij  IN  READ  ROUTINE 
c 

lu=80+ij-l 

write (lu) ( ( ( (rrij (j , i,iz,izp) , 

1  j=l, nyd2) , i=l,nx) , iz=l , nzd21) , izp=l,nzl) 


8000  CONTINUE 

write(6, 110) '  NORMALIZATION  FACTORS:' 

110  format (/a30) 

write (6, 111)  ( ( (ij , izp,cn(ij , izp) ) , izp=l,nzl) ,IJ=1,3) 

111  format (5 (2i6, flO. 6) ) 

stop 

end 

subroutine  setphi (a,b, IJ) 

*  EXPANDS  uiuj  OR  PHIij  ARRAYS  FROM  a(nym,nxm)  TO  b(ny,nxp2)  * 

*  i.  e.,  FROM  a(23,14)  TO  b(64,18)  * 

icitit'kititiei^-kie’kitic'k-kiiitieiiicie'k-kie'k’k'kit'kitit'k'k-k'k'kitic’k’kie’k’kit'k'k-k’k-k'kit'kit'k'k'k'k'k'kit'k'k'k’k'kicit’k'k'k'k’k 

parameter  (ny=64 , nx=16) 
parameter  (nym=23 , nxm=14) 
parameter  (nxp2=nx+2) 
dimension  a (nym, nxm) , b (ny , nxp2) 


do  1  i=l,nxm 

J 

1:14 

do  1  j=l,nym 

1 

1:23 

1 

b(j,i)=a(j/i) 

do  2  i=nxm+l,nxp2 

J 

15:18 

do  2  j=nym+i,33 

J 

24:33 

2 

b(j ,i)=0.0 

do  3  i=nxra+l,nxp2 

1 

15:18 

do  3  j=l,nym 

1 

1:23 

3 

o 

o 

II 

XI 

do  4  i=l,nxm 

J 

1,14 

do  4  j=nym+l,33 

1 

24:33 

4 

b(j,i)=0.0 

IF  (IJ.LE.4)  THEN 

1 

11,22, 

33,12 

do  10  i=l,nxp2 

1 

1:18 

do  10  j=34,ny 

J 

34:64 

10 

b(j,i)=  b(ny-j  +  2,i) 

ENDIF 

IF  (IJ.GE.5)  THEN 

J 

13,23 

have  odd  symmetry! 

do  11  i=l,nxp2 

J 

1:18 

do  11  j=34,ny 

J 

34:64 

11 

b(j , i)=-b(ny-j+2, i) 

ENDIF 
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return 

end 


subroutine  rij (a,b,cn, ij , iz, izp) 

*********************************************************************** 

*  INPUTS  EXPANDED  ARRAY  FROM  SETPHI,  OUTPUTS  CORRELATION  FUNCTN.  Rij  * 
*********************************************************************** 

parameter  (nz=32 , nzl=nz+l , ny=64 , nx=16) 
parameter  ( nxp2=nx+2 , nxd2 l=nx/2+l ) 
parameter  ( nyt 2=2 *ny , nyd2 l=ny/2+ 1 ) 

parameter  (nyd2=ny/2 , nxd2=nx/2,ny5=5*ny, nx3p4=3*nx+4) 
dimension  eyy(ny5) ,exx(nx3p4) 
dimension  vecy(nyt2) ,vecx(nxp2) 

dimension  a(ny,nxp2) ,b(nyd2 ,nxp2) , c ( nyd2 , nxp2 ) ,cn(3,nzl) 


c  INITIALIZE  FOURIER  TRANSFORMS  (CALC,  eyy  AND  exx) : 
c 

call  cfft2  (1, l,ny,vecy,eyy,vecy) 
call  crfft2 (1, l,nx,vecx,exx,vecx) 
c 

C  EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  Y  DIRECTION: 
c 

do  1  i=l,nxd21 
do  2  j=l,ny 

vecy(2*j-l)=a(j ,2*i-l)  !  REAL  PART 

2  vecy(2*j  )=a(j,2*i  )  !  IMAG  PART 

call  cfft2 (0, l,ny,vecy,eyy,vecy) 

do  3  j=l,nyd2  !  PUT  VECY  INTO  TEMPORARY  OUTPUT 

c(j , 2*i-l)=vecy (2*j-l) 

3  c(j,2*i  )=vecy(2*j  ) 

1  continue 
c 

C  EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  X  DIRECTION: 

C 

do  4  j=l,nyd2 

do  5  i=l,nxp2 

5  vecx(i)=c( j , i) 

call  crfft2 (0, l,nx, vecx,exx,vecx) 

do  6  i=l,nx  !  PUT  VECX  INTO  FINAL  OUTPUT 

6  c(j ,  i)=vecx(i) 

4  continue 

if  (IJ. LE. 3 . AND. izp. eq. iz)  then  !  set  normalization  factors 

cn(ij , iz)=c(l, 1)  !  ij=l,3  for  11,22,33 

if  (cn(ij , iz) .eq. 0. )  cn(ij,iz)=l. 
cn(ij ,nzl+i-iz)=cn(ij ,iz) 
endif 


c 
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C  FINALLY,  FLIP  RIGHT  SIDE  TO  THE  LEFT  OF  THE  LEFT  SIDE 
c 

do  10  j=l,nyd2 
do  11  i=l,nxd2 

b(j,i  )=c(j,i+8) 

11  b(j,i+8)=c(j,i  ) 

10  continue 

return 

end 
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program  wavspec 

************************************************************************ 

*  "Computes  space-time  corelations  for  planar  data"  * 

*  THIS  VERSION  READS  7  FILES,  DOES  5  FFT'S,  AND  AVOIDS  ZERO  PACKING  * 
************************************************************************ 

*  This  version  reads  three  3-d  yxt  files  ul,  u2,  and  u3,  extracts  from* 

*  each  three  2-d  arrays  (umn,unk,umk) ,  does  2-d  fft  on  each  of  those  * 

*  still  using  same  array  names,  and  adds  up  the  results  over  the  * 

*  remaining  (third)  variable.  It  then  calculates  energies,  ensemble  * 

*  averages,  and  computes  correlation  fundtions,  saving  all  * 

*  results  to  disk. 

parameter  (ny=64,nx=16,ntim=512) 

dimension  ul (ny , nx, ntim) ,u2 (ny,nx,ntim) ,u3 (ny,nx,ntim) 
dimension  umnl(ny  ,  nx+2) ,umn2 (ny  ,  nx+2) , umn3 (ny  ,  nx+2) 

dimension  smnl (ny/2 , nx/2+1) , smn2 (ny/2 , nx/2+1) , smn3 (ny/2 , nx/2+1) 
dimension  smnls(ny  ,  nx+2 ) , smn2s (ny  ,  nx+2) , smn3s (ny  ,  nx+2) 
dimension  rmnl(ny/2,  nx+2) , rmn2 (ny/2 ,  nx+2) , rmn3 (ny/2 ,  nx+2) 
dimension  unkl(nx  ,ntiin+2)  ,unk2  (nx  ,ntim+2)  ,unk3  (nx,  ntim+2) 
dimension  snkl (nx, ntim/2+1) ,snk2 (nx, ntim/2+1) , snk3 (nx, ntim/2+1) 
dimension  snkls(nx,  ntim+2) ,snk2s(nx,  ntim+2) , snk3s (nx,  ntim+2) 
dimension  rnkl (nx/2 , ntim+2) , rnk2 (nx/2 , ntim+2 ) ,rnk3 (nx/2 , ntim+2) 
dimension  umkl (ny, ntim+2)  ,umk2(ny,  ntim+2) ,umk3 (ny,  ntim+2) 
dimension  smkl (ny/2 , ntim/2+1) , smk2 (ny/2 , ntim/2+1) , 

*  smk3 (ny/2 , ntim/2+1) 

dimension  smkls (ny, ntim+2)  ,smk2s(ny, ntim+2)  ,smk3s(ny, ntim+2) 
dimension  rmkl (ny/2 , ntim+2) ,rmk2 (ny/2, ntim+2) ,rmk3 (ny/2, ntim+2) 
dimension  tmp(ny,nx) ,e(9) ,eavg(3) ,w(9) , tmpl (nx, ntim/2+1) 
c 

c  initialize  9  ensemble  averages  smni,  ski,  smki,  i=l,3  : 
c 

do  40  m=l,ny/2 
do  40  n=l, nx/2+1 
smnl (m, n) =0 . 
smn2 (m, n) =0 . 
smn3 (m, n) =0 . 

40  continue 

do  43  n=l,nx 
do  43  k=l, ntim/2+1 
snkl (n, k) =0 . 
snk2 (n,k)=0. 
snk3 (n,k) =0. 

43  continue 

do  46  m=l,ny/2 
do  46  k=l, ntim/2+1 
smkl (m, k) =0 . 
smk2 (m, k) =0 . 
smk3 (m,k)=0. 

46  continue 
c 

c  7  files  (lu=10,l6)  permit  5  FFT's  with  240  k's  left  over  at  the  end 


call  read  (ul , u2 , u3 , tmp, nx, ny ,ntim, 10 ,  1,400)  !  400 

call  read  (ul ,u2 , u3 , tmp, nx, ny , ntim, 11 , 401 , 512)  !  112 

call  hf ens  (ul , u2 , u3 , umnl , umn2 , umn3 , smnl , smn2 , smn3 , unkl , unk2 , unk3 , 
*  snkl , snk2 , snk3 , umkl , umk2 , umk3 , smkl , smk2 , smk3 , nx, ny , ntim)  1 ==== 

call  read  (ul,u2 ,u3 , tmp,nx,ny,ntim, 11,  1,288)  !  288 

call  read  (ul,u2,u3,tmp,nx,ny,ntim,12,289,512)  !  224 
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call  hf ens  (ul , u2 , u3 , umnl , umn2 , umnS , smnl , smn2 , smn3 , unkl , unk2 , unk3 , 

*  snkl ,  snk2 ,  snk3  ,  umkl ,  uink2 ,  umk3 ,  smkl ,  sink2  ,  sink3 ,  nx ,  ny ,  nt im)  !  ===== 

call  read  (ul, u2  ,u3 ,  tinp,nx,ny,ntim,  12 ,  1,176)  !  176 

call  read  (ul,u2,u3,tinp,nx,ny,ntiin,  13, 177, 512)  !  336 

call  hf ens  ( ul , u2 , u3 , umnl , umn2 , umn3 , smnl , smn2 , smn3 , unkl , unk2 , unk3 , 

*  snkl , snk2 , snk3 , umkl , umk2 , umk3 , smkl , smk2 , smk3 , nx , ny , ntim)  ! ==== 

call  read  (ul, u2 , u3 , tmp, nx,ny,ntim, 13 ,  1,  64)  !  64 

call  read  (ul,u2 ,u3 , tmp,nx,ny , ntim, 14 ,  65,464)  !  400 

call  read  (ul,u2 ,u3 , tmp,nx,ny ,ntim, 15 , 465 , 512)  !  48 

call  hf ens  (ul , u2 , u3 , umnl , umn2 , umn3 , smnl , smn2 , smn3 , unkl , unk2 , unk3 , 

*  snkl , snk2 , snk3 , umkl , umk2 , umk3 , smkl , smk2 , smk3 , nx , ny , ntim)  ! ==== 

call  read  (ul,u2,u3, tmp,nx,ny,ntim, 15,  1,352)  !  352 

call  read  (ul,u2 ,u3 , tmp,nx,ny,ntim, 16, 353 , 512)  !  160 

call  hf ens  (ul , u2 , u3 , umnl , umn2 , umn3 , smnl , smn2 , smn3 , unkl , unk2 , unk3 , 

*  snkl ,  snk2 ,  snk3  , umkl ,  umk2 , umk3 ,  smkl ,  smk2 ,  smk3 ,  nx,  ny ,  ntim)  !  ===== 

c 

c  Calculate  and  print  energies: 

e(l),w(l)) 
e(2) ,w(2)) 
e(3),w(3)) 

I 

!  average  over  5*512  data  pts 

I 

2,  3  are:  ' , (e(i) ,i=l,3) 

call  energ2(snkl,  nx,  ntim,  e(4),w(4)) 
call  energ2(snk2,  nx,  ntim,  e(5),w(5)) 
call  energ2(snk3,  nx,  ntim,  e(6),w(6)) 
e(4)=sqrt(e(4)/5./64.*512./512.)  ! 

e(5)=sqrt(e(5)/5./64.*512./512.)  I  energ2  divides  by  512**2 
e(6)=sqrt(e(6)/5./64.*512./512.)  1 

print  *,  '  NK  RMS  values  4,  5./,  6  are:  ' , (e(i) , i=4,6) 

call  energ3(smkl,  ny,  ntim,  e(7),w(7)) 
call  energ3(smk2,  ny,  ntim,  e(8),w(8)) 
call  energ3(smk3,  ny,  ntim,  e(9),w(9)) 
e(7)=sqrt(e(7)/5./16.*512./512. )  ! 

e(8)=sqrt(e(8)/5./16**512./512.)  !  energ3  divides  by  512**2 

e(9)=sqrt(e(9)/5./16.*512./512. )  ! 

print  *,  '  MK  RMS  values  7,  8,  9  are:  ' , (e (i) , i=7 , 9) 
c 

c  test  energies: 
c 

do  1050  i=l,3 
eavg(i)=0. 
do  1010  j=i,9,3 

eavg ( i ) =eavg ( i ) +e ( j ) 

1010  continue 

eavg ( i ) =eavg ( i ) / 3 . 
do  1020  j=i,9,3 

if (abs(e(j ) -eavg(i) )/eavg(i) .gt. 0.00001)  then 

print  *,*  energy  error!  e(',j,‘)  =  ',e(j),';  avg  =  •,eavg(i) 

end  if 

1020  continue 
1050  continue 
c 

c  normalize  ensemble  averages  by  energy: 
c 

do  1100  m=l,ny/2 


call  energl(smnl,  ny,  nx, 
call  energl(smn2,  ny,  nx, 
call  energl(smn3,  ny,  nx, 
e(l)=sqrt(e(l)/5./512.) 
e(2)=sqrt(e(2)/5./512. ) 
e(3)=sqrt(e(3)/5./512.) 
print  *,  '  MN  RMS  values  1, 
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do  1100  n=l,nx/2+l 

smnl  (m,  n)  =smnl  (in,n)/w(l) 
smn2  (in,n)=sinn2  (in,n)/w(2) 
smn3  (m,n)=sinn3  (in,n)/w(3) 

1100  continue 

do  1150  n=l,nx 
do  1150  k=l,ntim/2+l 

snkl ( n , k) =snkl ( n , k) /w ( 4 ) 
snk2 (n,k) =snk2 (n,k)/w(5) 
snk3 (n,k) =snk3 (n, k)/w(6) 

1150  continue 

do  1200  in=l,ny/2 
do  1200  k^l, ntim/2+1 

smkl(m,k)=sinkl(m,k) /w(7) 
sink2  (in,k)=sink2  (m,k)/w(8) 
sink3  (in,k)  =smk3  (ni,k)/w(9) 

1200  continue 
c 

c  symmetrize  smni,  snki,  smki: 
c 

call  setphil (smnl, smnls,ny,  nx) 
call  setphil (smn2 , smn2s,ny,  nx) 
call  setphil (smn3,smn3s,ny,  nx) 
call  setphi2 (snkl, snkls,nx,ntim) 
call  setphi2 (snk2 ,snk2s,nx,ntim) 
call  setphi2 (snk3,snk3s,nx,ntim) 
call  setphil (smkl, smkls,ny,ntim) 
call  setphil (smk2 , smk2s,ny,ntim) 
call  setphil (smk3,smk3s,ny,ntim) 
c 

c  calculate  correlation  functions: 
c 

call  rijl(sronls,rmnl,ny,  nx) 
call  rijl(smn2s,rmn2,ny,  nx) 
call  rij 1 (smn3s, rmn3 ,ny,  nx) 
call  rij2 (snkls,rnkl,nx,ntim) 
call  rij2 (snk2s, rnk2 ,nx,ntim) 
call  rij2 (snk3s, rnk3 ,nx,ntim) 
call  rij3 (smkls,rmkl,ny,ntim) 
call  rij3 (smk2s,rmk2,ny,ntim) 
call  rij3(smk3s,rmk3,ny,ntiro) 
c 

c  flip  snk  spectra  now  that  fft's  are  completed 
c 

do  1310  n=l,nx 
do  1310  k=l,ntim/2+l 
tmpl ( n , k) =snkl ( n , k) 

1310  continue 

do  1320  n=l,nx/2 
do  1320  k=l,ntim/2+l 

snkl(n,k)=tmpl(n+nx/2,k) 

1320  continue 

do  1330  n=nx/2+l,nx 
do  1330  k=l,ntim/2+l 

snkl (n,k)=tmpl (n-nx/2 ,k) 

1330  continue 

do  1410  n==l,nx 
do  1410  k=l,ntim/2+l 
tmpl (n,k)=snk2 (n,k) 

1410  continue 
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do  1420  n=l,nx/2 
do  1420  k=l,ntiin/2+l 

snk2 (n,k) =tmpl (n+nx/2,k) 

1420  continue 

do  1430  n=nx/2+l,nx 
do  1430  k=l,ntiin/2+l 

snk2  (n,  k)  =tinpl  (n-nx/2  ,k) 

1430  continue 

do  1510  n=l,nx 
do  1510  k=l,ntiin/2+l 
tmpl (n,k) =snk3 (n,k) 

1510  continue 

do  1520  n=l,nx/2 
do  1520  k=l,ntim/2+l 

snk3  (n,  k)  =tinpl  (n+nx/2  ,k) 

1520  continue 

do  1530  n=nx/2+l,nx 
do  1530  k=l,ntim/2+l 

snk3  (n,  k)  =tinpl  (n-nx/2  ,k) 

1530  continue 
c 

c  save  flipped  (where  applicable)  spectra  to  disk: 
c 

write(50,91)  l,e(l) ,eavg(l) 

write(50,92)  ( (smnl  (m,n)  ,  Tn=l,ny/2)  ,  n=l,  nx/2+1) 

write(50,91)  2,e(2) ,eavg(2) 

write(50,92)  ( (smn2  (m,n)  ,  in=l,ny/2)  ,  n=l,  nx/2+1) 

write(50,91)  3,e(3) ,eavg(3) 

write(50,92)  ((smn3(m,n),  in=l,ny/2),  n=l,  nx/2+1) 

write (50, 91)  4 ,e(4) ,eavg(l) 

write(50,92)  ((snkl(n,k),  n=l,  nx)  ,  k=l,ntiin/2+l) 
write(50,91)  5,e(5) ,eavg(2) 

write(50,92)  ((snk2(n,k),  n=l,  nx) ,  k=l,ntim/2+l) 
write(50,91)  6,e(6) ,eavg(3) 

write(50,92)  ((snk3(n,k),  n=l,  nx)  ,  k=l,ntiin/2+l) 

write(50,91)  7,e(7) ,eavg(l) 

write(50,92)  ( (smkl (m,k) ,  in=l,ny/2),  k=l, ntiin/2+1) 

write(50,91)  8,e(8) ,eavg(2) 

write (50, 92)  ( (smk2  (ra,  k)  ,  ni=l,ny/2)  ,  k=l,  ntiin/2+1) 

write(50,91)  9,e(9) ,eavg(3) 

write(50,92)  ( (sink3  (in,k)  ,  in=l,ny/2)  ,  k=l,ntiin/2+l) 

91  forinat(2x,i3,2el2.5) 

92  forTnat(2x,4el2.5) 
c 

c  save  correlation  functions; 
c 

write(50,91)  1 

write(50,92)  ( (rmnl (m, n) , m=l, ny/2 ) , n=l ,  nx) 

write(50,91)  2 

write(50,92)  ( (nnn2  (in,n)  ,in=l,ny/2)  ,n=l,  nx) 

write(50,91)  3 

write (50, 92)  ( (rmn3  (in,n)  ,in=l,ny/2)  ,n=l,  nx) 
write(50,91)  4 

write(50,92)  ( (rnkl(n,k)  ,n=l,nx/2)  ,k=l,ntini) 

write(50,91)  5 

write (50, 92)  vi^nk2  (n,k)  ,n=l,nx/2)  ,k=l,ntim) 
write(50,91)  6 

write (50, 92)  ( (rnk3  (n,k)  ,n=l,nx/2)  ,k=l,ntini) 

write(50,91)  7 
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write(50,92) 
write(50,91) 
write (50, 92) 
write (50,91) 
write(50,92) 


( (rmkl(m,k)  ,iii=l,ny/2)  ,k=l,ntim) 
8 

( (rink2  (m,k)  ,in=l,ny/2)  ,k=l,ntiin) 
9 

( (rmkO  (in,k)  ,in=l,ny/2)  ,k=l,ntiin) 


stop 

end 
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subroutine  read ( ul , u2 , u3 , tmp , nx , ny , nt im , lu , ka , kz ) 

dimension  ul (ny , nx, ntim) , u2 (ny , nx, ntim) ,u3 (ny,nx,ntim) ,tmp(ny,nx) 
do  50  k=ka,kz 

read  (lu)  ( (ul( j , i,k) , j=l,ny) , i=l,nx) 
read  ^u)  ( (u3  ( j  ,  i,k)  ,  j=l,ny) ,  i=l,nx) 

read  Uu)  Uu2  (j  ,  i,k)  ,  j=l,ny)  ,  i=l,nx) 
read  ^u)  ((tmp(j,i)  ,  j=l,  ny)  ,  i=l ,  nx) 
read  (lu)  ((tmp(j,i)  , j=l, ny) , i=l, nx) 
read  (lu)  ((tmp(j,i)  , j=l,ny) , i=l,nx) 
continue 
return 
end 


subroutine  hfens  (ul,u2,u3,  umnl,umn2 ,umn3 , smnl, smn2 , smn3 , 

*  unkl , unk2 , unk3 , snkl , snk2 , snk3 , umkl , umk2 , umk3 , smkl , smk2 , smk3 , 

*  nx,ny,ntim) 


dimension 

dimension 

dimension 

dimension 

dimension 

dimension 

dimension 


ul(ny,nx,ntim) ,u2 (ny,nx,ntim) ,u3 (ny,nx,ntim) 
umnl(ny  ,  nx+2) ,umn2 (ny  ,  nx+2) ,umn3 (ny  ,  nx+2) 

smnl (ny/2 , nx/2+1) , smn2 (ny/2 , nx/2+1) , smn3 (ny/2 , nx/2+1) 
unkl(nx  ,ntim+2) ,unk2 (nx  ,ntim+2) ,unk3 (nx,  ntim+2) 
snkl (nx, ntim/2+1) , snk2 (nx, ntim/2+1) , snk3 (nx, ntim/2+1) 
umkl (ny, ntim+2)  ,umk2(ny,  ntim+2) ,umk3 (ny,  ntim+2) 
smkl (ny/2 , ntim/2+1) ,smk2 (ny/2, ntim/2+1) , 
smk3 (ny/2 , ntim/2+1) 


c  At  each  time  (3rd  variable) ,  extract  umn  array,  fft  it,  compile  sum: 
c 

do  100  k=l,ntim 
do  70  m=l,ny 
do  70  n=l,nx 

umnl (m,n) =ul (m,n,k) 
umn2  (m,  n)  =u2  (m,  n, k) 
umn3  (m,  n)  =u3  (m,  n, k) 

70  continue 

call  hfftl (umnl,ny,nx, 1) 
call  ensavl (umnl, smnl, ny,nx) 
call  hfftl (umn2, ny, nx, 1) 
call  ensavl (umn2,smn2,ny,nx) 
call  hf ftl (umn3 , ny , nx, 1) 
call  ensavl (umn3,smn3,ny,nx) 

100  continue 
c 

c  At  each  y  (3rd  variable),  extract  unk  array,  fft  it,  compile  sum: 
c 

do  200  m=l,ny 
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do  170  n=l,nx 
do  170  k=l,ntiin 

unkl(n,k)=ul (m,n,k) 
unk2  (n,k)  =u2  (in,n,k) 
unk3(n,k)=u3  (in,n,k) 

170  continue 

call  hfft2  (unkl,nx,ntiin,  1) 
call  ensav2 (unkl,snkl,nx,ntim) 
call  hfft2 (unk2 , nx, ntim, 1) 
call  ensav2  (unk2 ,  snk2 ,  nx,ntiin) 
call  hfft2(unk3,nx,ntim,l) 
call  ensav2  (unk3,snk3,nx,ntiin) 

200  continue 
c 

c  At  each  x  (3rd  variable) ,  extract  umk  array,  fft  it,  compile  sum; 
c 

do  300  n=l,nx 
do  270  m=l,ny 
do  270  k=l,ntim 

umkl (m,k) =ul (m,n,k) 
umk2 (m,k)=u2 (m,n,k) 
umk3 (m,k)=u3 (m,n,k) 

270  continue 

call  hfft3 (umkl, ny, ntim, 1) 
call  ensav3 (umkl, smkl,ny, ntim) 
call  hfft3 (umk2 , ny,ntim, 1) 
call  ensav3 (umk2,smk2,ny,ntim) 
call  hfft3(umk3,ny,ntim,l) 
call  ensav3 (umk3,smk3,ny,ntim) 

300  continue 
return 
end 


subroutine  hfftl(a,nl,n2 , is) 

parameter  (ml=64 ,m2=16)  !  needed  to  dimension  noncalled  arrays 
dimension  a(nl,n2+2) 
dimension  indx(m2+2) , indy (2*ml) 
dimension  exx(3*m2+4) ,exy(5*ml) 
dimension  vecx(m2-t-2)  ,vecy(2*ml)  ,vecyl(2*ml) 
if  (is.eq.  -1)  go  to  500 
call  rcfft2 (1, l,n2 , vecx,exx,vecx) 
call  cfft2 (1, l,nl,vecy,exy,vecy) 
do  1  j=l,nl 
do  2  i=l,n2+2 

indx(i)=(i-l) *nl+j 
2  continue 

call  gather (n2 ,vecx, a, indx) 
call  ref ft2 (0,1, n2 , vecx, exx, vecx) 
call  scatter (n2+2, a, indx, vecx) 

1  continue 

do  3  i=l,n2/2+i 
do  4  j=l,nl 

jl=2*(j-l)+l 
indy{jl)=j+2*(i-l) *nl 
indy ( j 1+ 1 ) =indy ( j 1 ) +nl 
4  continue 

call  gather (2*nl, vecy, a, indy) 
call  cfft2 (0, l,nl,vecy,exy,vecy) 
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call  scatter(2*nl,a, indy,vecy) 

3  continue 
return 

500  continue 

call  erf ft2 (1, -1, n2 , vecx,exx, vecx) 
call  cfft2 (l,-l,nl,vecy,exy,vecy) 
do  5  i=l,n2/2+l 
do  6  j=l,nl 

jl=2*(j-l)+l 
indy ( jl)=j+2* (i-1) *nl 
indy ( j 1+1) =indy ( j 1) +nl 

6  continue 

call  gather (2*nl , vecy , a, indy) 
call  cfft2 (0, -1, nl, vecy , exy, vecy) 
call  scatter ( 2 *nl, a, indy, vecy) 

5  continue 

do  7  j=l,nl 
do  8  i=l,n2+2 

indx(i)=(i-l) *nl+j 
8  continue 

call  gather (n2+2 , vecx, a, indx) 
call  erf ft2 (0 , -1, n2 , vecx, exx, vecx) 
call  scatter (n2 , a, indx, vecx) 

7  continue 

do  10  i=l,n2 
do  10  j=l,nl 

ap  ,  i)=a(  j  ,  i)/  (2  .  *float  (n2*nl)  ) 

10  continue 

return 
end 


subroutine  hf ft2 (a, nl , n2 , is) 

parameter  (ml=16,m2=512)  !  needed  to  dimension  noncalled  arrays 
dimension  a(nl,n2+2) 
dimension  indx(m2+2) , indy (2*ml) 
dimension  exx(3*m2+4) ,exy(5*ml) 
dimension  vecx(m2+2) ,vecy(2*ml) ,vecyl(2*ml) 
if  (is. eg.  -1)  go  to  500 
call  rcfft2 (1, l,n2 , vecx, exx, vecx) 
call  cfft2 (1, l,nl, vecy, exy, vecy) 
do  1  j=l,nl 
do  2  i=l,n2+2 

indx (i)=( i-1) *nl+j 
2  continue 

call  gather (n2 , vecx, a, indx) 
call  rcfft2 (0, l,n2, vecx, exx, vecx) 
call  scatter (n2+2 , a, indx, vecx) 

1  continue 

do  3  i=l,n2/2+l 
do  4  j=l,nl 

jl=2*(j-l)+l 
indy ( j 1) =j+2* (i-1) *nl 
indy ( j 1+1 ) =indy ( j 1 ) +nl 
4  continue 

call  gather (2 *nl, vecy, a, indy) 
call  cfft2 (0, l,nl, vecy, exy, vecy) 
call  scatter(2*nl,a, indy,vecy) 

3  continue 
return 
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500  continue 

call  crfft2 (1, -l,n2,vecx,exx,vecx) 
call  cfft2 (1, -l,nl,vecy,exy,vecy) 
do  5  i=l,n2/2+i 
do  6  j=l,nl 

jl=2*(j-l)+l 
indy ( j 1) =j+2* ( i-1) *nl 
indy (j 1+1) =indy(jl)+nl 

6  continue 

call  gather ( 2*nl , vecy , a , indy) 
call  cfft2 (0,-l,nl, vecy, exy, vecy) 
call  scatter (2 *nl, a, indy, vecy) 

5  continue 

do  7  j=l,nl 
do  8  i=l,n2+2 

indx ( i) = ( i-1) *nl+j 
8  continue 

call  gather (n2+2 ,vecx, a, indx) 
call  crfft2 (0 , -1 , n2 , vecx,exx, vecx) 
call  scatter (n2 , a, indx, vecx) 

7  continue 

do  10  i=l,n2 
do  10  j=l,nl 

a(j,i)=a(j,i)/(2. * float (n2*nl) ) 

10  continue 

return 
end 


subroutine  hfft3 (a,nl,n2, is) 

parameter  (rol=64 ,m2*512)  !  needed  to  dimension  noncalled  arrays 
dimension  a(nl,n2+2) 
dimension  indx (m2+2) , indy (2*ml) 
dimension  exx(3*m2+4) ,exy (5*ml) 
dimension  vecx(m2+2) , vecy (2*ml) , vecyl (2*ml) 
if  (is.eq.  -1)  go  to  500 
call  ref ft2 (1 , 1, n2 , vecx,exx,vecx) 
call  cfft2 (1, l,nl, vecy, exy, vecy) 
do  1  j=l,nl 
do  2  i=l,n2+2 

indx ( i) = (i-1) *nl+j 
2  continue 

call  gather (n2 , vecx, a, indx) 
call  rcfft2 (0, l,n2,vecx,exx,vecx) 
call  scatter (n2+2 , a, indx, vecx) 

1  continue 

do  3  i=l,n2/2+l 
do  4  j=l,nl 

jl=2*(j-l)+l 
indy ( j 1) =j+2* ( i-1) *nl 
indy (j 1+1) =indy(jl)+nl 
4  continue 

call  gather (2 *nl, vecy, a, indy) 
call  cfft2 (0, l,nl, vecy, exy, vecy) 
call  scatter(2*nl,a, indy,vecy) 

3  continue 
return 

500  continue 

call  crfft2 (1, -l,n2, vecx,exx,vecx) 
call  cfft2 (l,-l,nl, vecy, exy, vecy) 
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do  5  i=l,n2/2+l 
do  6  j=l,nl 

jl=2*(j-l)+l 
indy( jl)=j+2* (i-1) *nl 
indy ( j 1+1) =indy ( j 1) +nl 
continue 

call  gather (2*nl,vecy,a, indy) 
call  cfft2 (0, -l,nl, vecy,exy,vecy) 
call  scatter (2*nl , a, indy jVecy) 

5  continue 

do  7  j=l,nl 
do  8  i=l,n2+2 

indx ( i) = (i-1) *nl+j 
8  continue 

call  gather (n2+2,vecx, a, indx) 
call  crfft2 (0 , -1 , n2 , vecx,exx, vecx) 
call  scatter (n2 , a, indx, vecx) 

7  continue 

do  10  i=l,n2 
do  10  j=l,nl 

a(j , i)=a(j , i) / (2 . *float (n2*nl) ) 

10  continue 

return 
end 


subroutine  ensavl (a,b,nl,n2) 
parameter  (ml=64 ,m2=16) 
dimension  a(nl,n2+2) ,b(nl/2,n2/2+l) 

dimension  r(ml,m2/2+l)  !  need  var  dim  =  r(nl,n2/2+l) 
do  1  i=l,n2/2+l 

do  1  j=l,nl 

r(j,i)=a(j,2*i-l) **2+a ( j , 2*i) **2 

1  continue 

call  syml (r,nl,n2/2+l) 
do  2  i=l,n2/2+l 
do  2  j=l,nl/2 

b(j,i)=b(j,i)+r(j,i) 

2  continue 
return 
end 


subroutine  ensav2 (a,b,nl,n2) 
dimension  a(nl,n2+2) ,b(nl,n2/2+l) 
do  2  i=l,n2/2+l 

do  2  j=l,nl 

b(j,i)=b(j,i)+a(j,2*i-l) **2+a ( j , 2*i) **2 
2  continue 
return 
end 


subroutine  ensavl (a, b,nl,n2) 
parameter  (ml=64 ,m2=512) 
dimension  a(nl,n2+2) ,b(nl/2,n2/2+l) 

dimension  r(ml,m2/2+l)  !  need  var  dim  =  r(nl,n2/2+l) 
do  1  i=l,n2/2+l 
do  1  j=l,nl 

r(j , i)=a(j ,2*i-l) **2+a ( j , 2*i) **2 
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continue 

call  syml (r , nl , n2/2+l) 
do  2  i=l,n2/2+l 
do  2  j=l,nl/2 

b(j/i)=t'(j/i)+r(j,i) 

2  continue 

return 
end 


subroutine  syml (c, ml, m2) 
dimension  c(ml,m2) 
do  1  i=l,m2 
do  1  j=l,ml/2 

if  (  j  .eq.  1)  then 
c(j,i)  =  c(j,i) 
else 

c(j,i)  =  (c(j,i)  +  c(ml-j+2, i) )/2. 
end  if 

1  continue 
return 
end 


subroutine  energl(b,nl,n2,e,w) 

parameter  (ml=64 ,m2=16)  !  needed  to  dimension  noncalled  arrays 
dimension  b (nl/2 , n2/2+i) 
dimension  c (ml/2 ,m2/2+l) 
do  4  j=2,nl/2 

c(j,l)=b(j,l)/2. 

b(j,l)=b(j,l)/8. 

4  continue 

do  5  i=2,n2/2+l 
c(l,i)=b(l,i)/2. 
b(l,i)=b(l,i)/8. 

5  continue 
c(l, 1)=0. 
b(l,l)=b(l,l)/16. 
do  6  i=2,n2/2+l 
do  6  j=2,nl/2 

c(j,i)=b(j,i) 

b(j,i)=b(j,i)/4. 

6  continue 
e=0. 
w=0. 

do  7  i=l,n2/2+l 
do  7  j=l,nl/2 
w=w+c( j , i) 
e=e+b(j,i) 

7  continue 

e=e*4 ./nl/nl/n2/n2 

return 

end 


subroutine  energ2 (b,nl,n2,e,w) 

parameter  (ml=16,m2»512)  !  needed  to  dimension  noncalled  arrays 
dimension  b(nl,n2/2+l) 
dimension  c(ml,m2/2+l) 
do  4  j=2,nl 
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c(j,l)=b(j,l)/4. 

b(j,l)=b(j,l)/16. 

4  continue 

do  5  i=2,n2/2+l 
c(l,i)=b(l,i)/2. 
b(l,i)=b(l,i)/8. 

5  continue 
c(l,l)=0. 
b(l,l)=b(l,l)/16. 
do  6  i=2,n2/2+l 
do  6  j=2,nl 

c(  j  ,i)=b(j  ,i)/2. 
b(j  ,i)=b(j  ,  i)/8. 

6  continue 
w=0. 
e=0. 

do  7  i=l,n2/2+l 
do  7  j=l,nl 
w=w+c(j,i) 
e=e+b(j , i) 

7  continue 

e=e*4 ./nl/nl/n2/n2 

return 

end 


subroutine  energS (b,nl,n2,e,w) 

parameter  (ml=64 ,m2=512)  i  needed  to  dimension  noncalled  arrays 
dimension  b(nl/2,n2/2+l) 
dimension  c (ml/2, m2/2+l) 
do  4  j=2,nl/2 

c(j,l)*b(j,l)/2. 

b(j,l)=b(j,l)/8. 

4  continue 

do  5  i=2,n2/2+l 
c(l,i)=b(l,i)/2. 
b(l,i)=b(l,i)/8. 

5  continue 
c(l,l)=0. 
b(l,l)=b(l,l)/16, 
do  6  i=2,n2/2+l 
do  6  j=2,nl/2 

c(j,i)=b(j,i) 

b(j,i)=b(j,i)/4. 

6  continue 
e=0 . 
w=0. 

do  7  i=l,n2/2+l 
do  7  j=l,nl/2 
w=w+c ( j , i ) 
e=e+b(j,i) 

7  continue 

e=e*4 ./nl/nl/n2/n2 

return 

end 


subroutine  setphil (a,b,nl,n2) 
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o  o 


*  EXPANDS 

smni 

ARRAYS 

FROM 

a(ny/2,nx/2+l) 

TO 

b(ny,nx+2) 

* 

* 

i.  e. , 

FROM 

a(32. 

9) 

TO 

b(64,  18) 

* 

* 

smki 

ARRAYS 

FROM 

a (ny/2 , ntim/2+1) 

TO 

b(ny, ntim+2) 

* 

* 

i .  e.  , 

FROM 

a(32. 

257) 

TO 

b(64,514) 

* 

h*********************************************************************** 

dimension  a (nl/2 , n2/2+l) ,b(nl,n2+2) 

do  1  i=l,n2/2+l  !  1:  9,257 

do  1  j=l,nl/2  !  1:  32,  32 

b(j,2*i-l)  =a(j,i)  !  a  into  real(b) 
b(j,2*i  )  =0.  !  imag(b)=0. 

1  continue 

do  2  i=l,n2+2  1  1:  18,514 

b(nl/2+l, i) =0.  !  row  33,  33  =  0. 

2  continue 

do  10  i=l,n2+2  !  1:  10,  514 

do  10  j=nl/2+2,nl  !  34:64,  34:64 

b(j,i)=  b(nl-j+2,i)l  symmetry  for  smni,  smki 
10  continue 

b(l,l)=0.  !  elimininate  (0,0)  mode 

return 

end 


subroutine  setphi2 (a,b,nl,n2) 

*  EXPANDS  snki  ARRAYS  FROM  a(nx,  ntim/2+1)  TO  b(nx,  ntim+2)  * 

*  i.  e.,  FROM  a(16,  257)  TO  b(16,  514  )  * 

dimension  a(nl,n2/2+l) ,b(nl,n2+2) 

do  1  i=l,n2/2+l  I  1:  257 

do  1  j=l,nl  !  1:  16 

b(j,2*i-l)  =a(j,i)  !  a  into  real(b) 

b(j,2*i  )  =0.  !  imag(b)=0. 

1  continue 

b(l,l)=0.  !  elimininate  (0,0)  mode 

return 

end 


subroutine  rijl(a,b,nl,n2) 

************************************11********************************** 

*  INPUTS  EXPANDED  ARRAY  FROM  SETPHl,  OUTPUTS  CORRELATION  FUNCTN.  Rij  * 
*********************************************************************** 

parameter  (ml=64 ,m2=16)  !  needed  for  uncalled  array  dimensions 

dimension  eyy (5*ml) ,exx(3*m2+4) 
dimension  vecy (2*ml) , vecx(m2+2) 
dimension  a(nl,n2+2) ,b(nl/2,n2+2) ,c(ml/2 ,m2+2) 

INITIALIZE  FOURIER  TRANSFORMS  (CALC,  eyy  AND  exx) : 
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call  cfft2  (1, l,nl,vecy,eyy,vecy) 
call  crfft2 (1, l,n2,vecx,exx,vecx) 


EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  Y  DIRECTION: 
do  1  i=l,n2/2+l 
do  2  j=l,nl 

vecy(2*j-l)=a(j ,2*i-l)  !  REAL  PART 

2  vecy(2*j  )=a(j,2*i  )  !  IMAG  PART 

call  cfft2 (0, l,nl,vecy,eyy,vecy) 


do  3  j=l,nl/2 

c(j ,2*i-l)=vecy(2*j-l) 
3  c(j,2*i  )=vecy(2*j  ) 

1  continue 


!  PUT  VECY  INTO  TEMPORARY  OUTPUT 


EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  X  DIRECTION; 


do  4  j=l,nl/2 

do  5  i=l,n2+2 
5  vecx(i) =c( j , i) 

call  crfft2 (0, l,n2 ,vecx,exx,vecx) 


do  6  i=l,n2 
6  c( j , i) =vecx(i) 

4  continue 


!  PUT  VECX  INTO  FINAL  OUTPUT 


cnorm=c(l, 1) 

if  (cnorm.eq.  0. )  cnorTn=l. 

FINALLY,  FLIP  RIGHT  SIDE  TO  THE  LEFT  OF  THE  LEFT  SIDE; 

do  10  j=l,nl/2 
do  11  i=l,n2/2 

b(j,i  )=c(j  ,  i+n2/2)/cnorin 

11  b( j , i+n2/2) =c( j , i  )/cnonn 

10  continue 

return 

end 


subroutine  rij2 (a,b,nl,n2) 

1c********************************************************************** 

*  INPUTS  EXPANDED  ARRAY  FROM  SETPHI,  OUTPUTS  CORRELATION  FUNCTN.  Rij  * 
*********************************************************************** 


parameter  (ml=16,m2=512)  !  needed  for  uncalled  array  dimensions 

dimension  eyy (5*ml) ,exx(3*m2+4) 

dimension  vecy(2*ml) ,vecx(m2+2) 

dimension  a(nl,n2+2) ,b(nl/2,n2+2) ,c(ml/2,m2+2) 

INITIALIZE  FOURIER  TRANSFORMS  (CALC,  eyy  AND  exx) : 
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call  cfft2  (1, l,nl,vecy,eyy,vecy) 
call  crfft2 (1, l,n2 ,vecx,exx,vecx) 
c 

c  EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  Y  DIRECTION: 
c 

do  1  i=l,n2/2+l 
do  2  j=l,nl 

vecy(2*j-l)=a(j,2*i-l)  !  REAL  PART 

2  vecy(2*j  )=a(j,2*i  )  !  IMAG  PART 

call  cfft2 (0, l,nl,vecy,eyy,vecy) 

do  3  j=l,nl/2  1  PUT  VECY  INTO  TEMPORARY  OUTPUT 

c(j ,2*i-l)=vecy (2*j-l) 

3  c(j,2*i  )=vecy(2*j  ) 

1  continue 
c 

c  EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  X  DIRECTION: 
c 

do  4  j=l,nl/2 

do  5  i=l,n2+2 

5  vecx(i)=c( j  ,  i) 

call  crfft2 (0, l,n2 ,vecx,exx,vecx) 

do  6  i=l,n2  !  PUT  VECX  INTO  FINAL  OUTPUT 

6  c ( j , i ) -vecx ( i ) 

4  continue 
cnorin=c  (1,1) 

if  (cnorm.eq.O.)  cnorm=l. 
c 

c  FINALLY,  FLIP  RIGHT  SIDE  TO  THE  LEFT  OF  THE  LEFT  SIDE: 
c 

do  10  j=l,nl/2 
do  11  i=l,n2/2 

b(j,i  )=c( j , i+n2/2)/cnorm 

11  b(j ,i+n2/2)=c(j , i  )/cnorm 
10  continue 

return 

end 


subroutine  rij3 (a,b,nl,n2) 

*********************************************************************** 
*  INPUTS  EXPANDED  ARRAY  FROM  SETPHI,  OUTPUTS  CORRELATION  FUNCTN.  Rij  * 
*****************************************************************11***** 

parameter  (ral=64 ,m2=512)  J  needed  for  uncalled  array  dimensions 
dimension  eyy (5*ml) ,exx(3*m2+4) 
dimension  vecy (2*ml) , vecx(ro2+2) 
dimension  a(nl,n2+2) ,b(nl/2,n2+2) ,c(ml/2 ,m2+2) 
c 

c  INITIALIZE  FOURIER  TRANSFORMS  (CALC,  eyy  AND  exx) ; 
c 

call  cfft2  (1 , l,nl, vecy , eyy, vecy) 


44 


call  crfft2 (1, l,n2,vecx,exx,vecx) 


c  EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  Y  DIRECTION; 
c 

do  1  i=l,n2/2+l 
do  2  j=l,nl 

vecy(2*j-l)=a(j ,2*i-l)  I  REAL  PART 

2  vecy(2*j  )=a(j,2*i  )  i  IMAG  PART 

call  cf ft2 (0, l,nl, vecy,eyy,vecy) 

do  3  j=l,nl/2  !  PUT  VECY  INTO  TEMPORARY  OUTPUT 

c( j , 2*i-l) =vecy (2*j-l) 

3  c(j,2*i  )=vecy(2*j  ) 

1  continue 
c 

c  EXTRACT  VECTOR  FROM  ARRAY  AND  DO  FFT  IN  X  DIRECTION: 
c 

do  4  j=l,nl/2 

do  5  i=l,n2+2 

5  vecx(i)=c(j  ,  i) 

call  crfft2 (0, l,n2 , vecx,exx,vecx) 

do  6  i=l,n2  !  PUT  VECX  INTO  FINAL  OUTPUT 

6  c(j , i)=vecx(i) 

4  continue 
cnonn=c  (1,1) 

if  (cnorm.eq. 0. )  cnonn=l. 


c  FINALLY,  FLIP  RIGHT  SIDE  TO  THE  LEFT  OF  THE  LEFT  SIDE: 
c 

do  10  j=l,nl/2 
do  11  i=l,n2/2 

b(j,i  )=c(j  ,i+n2/2)/cnorin 
11  b(j , i+n2/2)=c( j , i  )/cnonn 
10  continue 

return 

end 
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